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ABSTRACT 


A  numerical  procedure  for  the  structural  analysis  of  a  general  three- 
dimensional  nature  has  been  developed  to  provide  a  reliable  solution  to  the 
problem  of  determining  the  strength  of  propellers,  particularly  those  with 
unconventional  configurations.  A  finite  element  displacement  model  is  utilized 
and  compatible  solid  elements  in  their  general  form  are  adopted.  The  use  of 
interpolation  functions  to  define  pertinent  curvilinear  coordinates  in  element 
space  gives  the  finite  element  technique,  new  capabilities  for  dealing  with 
structures  of  highly  complex  geometry.  This  formulation  bypasses  the  con¬ 
straints  of  simplifying  assumptions  (such  as  those  imposed  by  classicial  plate 
theory)  and  allows  a  closer  approximation  to  the  true  structural  configuration 
than  is  possible  by  other  approaches,  including  most  analytical  and  numerical 
methods.  The  performance  of  the  refined  elements  described  in  this  report  is 
distinctly  superior  to  those  obtainable  with  commonly  available  elements,  for 
example,  those  in  NASTRAN.  A  highly  skewed  propeller  blade  under  pre¬ 
scribed  pressure  distributions  was  chosen  for  demonstration  of  the  generality 
of  the  procedure.  Good  agreement  was  obtained  with  measured  displacement 
and  experimental  stress  data. 

ADMINISTRATION  INFORMATION 
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CHAPTER  1 


INTRODUCTION 


1.1  General 

Recent  years  have  witnessed  an  increased  interest  in  the  development  of  performance- 
oriented  surface  ships  for  which  it  is  vital  to  keep  weight  to  a  minimum,  e.g.,  high-speed 
hydrofoils,  catamarans,  and  surface  effect  ships.  A  more  accurate  method  of  analysis  than 
currently  employed  in  shipyards  is  imperative  if  weight  saving  is  to  be  achieved  for  such 
vehicles.  The  effective  use  of  materials  and  an  increased  reliability  of  design  will  have 
far-reaching  results  over  their  life  spans. 

Governed  by  functional  requirements  and  hydrodynamical  considerations,  the  geometry 
of  ship  scantlings  is  generally  complex  and  their  construction  contains  a  high  degree  of 
redundancy.  It  is  therefore  necessary  to  make  simplifying  assumptions  in  order  to  reduce 
the  complexity  of  the  mathematical  model  representing  the  structure  to  a  form  that  is  amenable 
to  traditional  design  methods.  In  consequence,  certain  characteristic  behavior  of  the  elastic 
body  is  ignored  and  the  accuracy  of  the  analysis  is  often  open  to  question.  Thus  verification 
requires  testing  scaled  models  and,  at  times,  costly  prototypes  as  well. 

The  rapid  advances  in  digital  electronic  computers  since  the  mid-1950’s  coupled  with 
recent  development  in  discrete  element  methods  provide  a  powerful  new  tool  for  structural 
analysis  (Paulling,  1964;  Moe  and  Tonnesen,  1966;  and  the  International  Ship  Structures 
Committee,  1969).*  Many  complicated  design  problems  that  were  considered  insurmountable 
to  a  realistic  analysis  only  a  few  years  ago  can  now  be  executed  almost  routinely  by  using 
an  ordinary  computer  (Roren,  1969;  Ma,  1969;  and  Abrahamsen,  1970).  Specifically  a  struc¬ 
ture  system  having,  say,  1000  degrees  of  freedom  can  be  solved  in  a  matter  of  a  few  min¬ 
utes  on  a  late  model  computer  (such  as  CDC  6600,  or  IBM  360/75,  etc.)  through  appropriate 
idealization  with  due  consideration  for  the  bandwidth  of  the  resulting  system  of  equations. 

During  the  past  decade,  the  development  of  finite  element  methods  has  exhibited  an 
exponential  growth,  and  the  demand  for  appropriate  programs  has  increased  rapidly.  The 
total  number  of  finite  element  computer  programs  in  which  substantial  efforts  have  been  ex¬ 
pended  may  have  exceeded  several  hundred  (Gallagher,  1970  and  Schrem,  1971).  Neverthe¬ 
less,  only  a  rather  small  number  of  them  (and  these  only  recently)  are  accessible  to  engi¬ 
neers  in  practice.  Among  these  well  known  programs  are  NASTRAN,  the  NASA  structural 
analysis  (MacNeal  and  McCormick,  1967);  STRUDL,  structural  design  language  (Logcher 
and  Sturman,  1966);  FINEL  (Adamchak,  1970)  and  SAMIS  (Melosh  et  al. ,  1966).  Programs  of 
a  proprietary  nature  includes  ASKA  (Schrem  and  Roy,  1971);  DAISY  (Kamel  et  al.,  1969); 
SESAM-69  (Araldsen  and  Egeland,  1971);  SAP  (Wilson,  1970);  STARDYNE  (Dainora,  1971); 
and  others  (Hartung,  1970;  and  Mallett  and  Jordan,  1969). 

♦References  are  listed  alphabetically  starting  on  page  95. 
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Most  of  the  finite  element  programs,  including  larger  scale  and  general  purpose  pro¬ 
grams  that,  are  currently  available,  have  little  or  no  data-generation  features  and  their  ele¬ 
ment  libraries  contain  basically  one-  and  two-dimensional  elements  of  linear  or  constant- 
strain  type.s.  The  big  cost  item  of  a  finite  element  analysis  is  frequently  the  data  preparation 
stage.  In  areas  of  steep  stress  gradient,  for  example,  the  good  approximation  of  an  important 
structural  response  may  require  the  assemblage  of  a  large  number  of  elements,  especially 
when  the  elements  to  be  used  are  of  lower  calibre,  such  as  the  constant-stress  elements.  In 
the  case  of  an  irregular  boundary,  curved  elements  will  have  a  distinct  advantage.  Thus  re¬ 
finement  of  elem  ent  characteristics  can  have  a  profound  impact  on  the  economy  and  range  of 
solution  that  the  finite  element  method  can  provide. 

By  nature  of  their  geometric  proportions,  many  structural  components  can  be  idealized 
as  one-  or  two-dimensional  problems  of  elasticity  and  standard  methods  can  be  used  to  obtain 
reasonably  good  solutions.  For  other  design  tasks,  however,  no  conventional  approach  can 
achieve  realistic  re'sults.  Examples  are  bodies  of  complex,  unsymmetrical  shapes  and  inter¬ 
face  problems  of  two  or  more  geometrical  entitles  (e.g.,  the  junction  of  pipes,  plates,  and/or 
shells).  The  rational  solution  of  such  problems  requires  a  general  analysis  in  three- 
dimensional  elasticity.  It  is  in  this  difficult,  but  important,  field  of  three-dimensional  prob¬ 
lems  that  recent  developments  in  the  isoparametric  element  family  offer  the  most  promising 
approach.  These  refined  elements  will  be  utilized  in  the  study  reported  here. 

1.2  Objective  and  Scope 

The  objective  of  the  present,  study  was  to  develop  a  numerical  procedure  for  the  static 
analysis  of  a  three-dimensional  elastic  body  of  arbitrary  configuration.  More  specifically, 
the  purpose  was  to  determine  the  structural  behavior  of  a  marine  propeller  subjected  to  a  pre¬ 
scribed  pressure  loading.  A  highly  skewed  propeller*  was  chosen  to  demonstrate  the  gener¬ 
ality  of  the  approach.  The  study  employed  the  finite  element  method  in  conjunction  with 
curved  solid  elements,.  A  computer  program  was  developed  to  implement  the  procedure-  for 
predicting  displacements  and  stresses  of  a  complex  structure  with  reference  to  an  arbitrary 
curvilinear  coordinate  system.  Linear  elasticity  and  small  deformation  theory  were  assumed. 

Chapter  2  outlines  the  finite  element  method  for  structural  analysis.  Certain  element 
characteristics  derived  f.-om  a  displacement  model  are  discussed  to  aid  in  the  selection  of 
appropriate  elements  for  improved  computational  results. 

Although  finite  element  techniques  are  widely  used  in  the  two-dimensional  domain  of 
plates  and  shells,  they  h.ave  had  only  limited  application  for  the  treatment  of  complex  struc¬ 
tures  in  the  context  of  three-dimensional  elasticity.  The  principal  reason  for  this  slow 

*The  marine  propeller  takes  on  complex,  skewed  geometry  as  a  result  of  design  considerations,  such  as  those 
of  vibration  and  cavitation  aspects  in  blade  design  (Cox  and  Morgan,  1972). 
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progress  is  the  large  amount  of  input  data  and  processing  time  required  to  implement  a  three- 
dimensional  solution  when  only  simple,  tetrahedron-type  solid  elements  are  utilized.  Isopara¬ 
metric  formulation  and  its  associated  refined  curved  elements  coupled  with  a  more  efficient 
solution  technique  (as  described  in  Chapter  3)  now  make  it  possible  to  tackle  some  of  the 
most  difficult  problems  in  solid  mechanics. 

Some  selected  problems,  including  beams,  plates,  shells,  and  stiffened  plates  are 
solved  to  evaluate  the  adequacy  and  performance  of  the  procedure  developed  here.  Further, 
a  ship  component  of  complex  geometry— a  skewed  propeller  blade— is  analyzed  in  Chapter  4 
to  provide  insight  into  the  potential  of  the  procedure  as  a  design  tool.  Since  no  analytic 
solution  for  the  propeller  blade  problem  is  known,  computed  results  for  displacements  and 
stresses  are  compared  with  experimental  data. 

1.3  Notations 

The  symbols  used  in  this  study  are  defined  where  they  first  appear.  For  convenience, 
frequently  used  symbols  are  summarized  below. 

The  bar  and  tilde  underscores  (or  overscores)  generally  denote  a  vector  and  a  matrix, 
respectively.  Parentheses  and  brackets  are  used  alternatively  to  denote  a  vector  and  a 
matrix.  For  example,  a  column  vector  Va  can  be  written 


with  its  subvector 

(u. ' 

K 

4  matrix 

6  =  [0]  =  [Vv  V2,  V3) 

When  vector  quantities  appear  in  an  equation,  standard  vector  notations  will  apply.  For 
example,  “x”  and  represent  cross  and  scalar  product  of  vectors,  respectively. 
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[A] 

[B] 

[B'] 

D 

ID] 

W'l 

DK 

E 

\F\ 


matrix  of  functions  in  nodal  coordinates,  Eq.  (2.2) 

matrix  relating  strain  vector  to  nodal  displace¬ 
ments,  Eq.  (2.4)  or  (3.8) 

matrix  relating  local  strain  vector  of  a  shell  to 
nodal  displacements,  Eq.  (3.52) 

a  specified  domain,  such  as  a  given  volume  or  area 

elasticity  matrix 

elasticity  matrix  in  the  local  coordinate  system  for 
a  isotropic  shell  element 

E 

(1  +  v)  (1  -  2v) 

Young’s  modulus  of  elasticity 

equivalent  load  vector,  also  known  as  generalized 
load  vector 


Fx(f),Fy(I),Fz(I) 

G 


equivalent  element  nodal  force  in  the  x,  y,  or  z 
direction,  respectively;  forces  are  positive  in 
the  positive  direction  of  x ,  y,  and  z  axes 

equivalent  force  in  direction  of  x,  y  or  z  axis  at 
node  “i,”  Eq.  (3.28) 

matrix  containing  functions  of  curvilinear  coordinates 


E 

- -  shear  modulus  of  elasticity 

2(1  +  v) 


HiX >  Hir  Hi: 


weighting  coefficient  corresponding  to  position  along 
Gaussian  quadrature  points  £.x,  rj.  ,  or 


subscript  indicating  nodal  number,  or  active  index 


I 


moment  of  inertia  of  a  transverse  section  of  a  beam 


AAA 

i,  j,  k 


vector  having  unit  value  in  direction  of  x,  y,  or 
z  axis,  respectively 


[J] 


Jacobian  matrix  of  coordinate  transformation 


w\ 


Jacobian  determinant 


6 


k 


constant  factor  included  in  [ D ']  matrix  to  improve 
shear  deformation 


k.. 

stiffness 

m 

stiffness 

stiffness 

coefficient  at  i  th  row  and  j  th  column 

matrix  of  entire  structure 
matrix  of  an  element  e 


K  (t  ■  s) 


submatrices  of  [Ke\ 


_  AAA 

l,m,n  direction  cosines  of  a  unit  vector  e  and  (i,  j,  k) 

representing  global  rectangular  axes  (x,  y,  z) 


[La] 


localizing  matrix  relating  element  nodal  parameter 
to  global  structure  parameter,  Eq.  (2.13) 


M  M 


bending  moment  components  in  the  x  and  y  directions 


n 


vector  normal  to  a  curved  surface 


N .  (£,  rj),  N  .(f,  77,  C)  function  of  curvilinear  coordinates  in  two  or  three 

dimensions,  respectively,  taking  a  value  of  unity 
at  node  i  and  zero  at  all  other  nodes 

NNPE  number  of  nodes  per  element 


P 

(Pi 
\Pa  I 

pe ■  *V 

q(x,  y>  z) 
W 

r 

r. 

1 

[R] 

S 


applied  pressure  on  an  element  face,  Eq.  (3.27) 
global  load  vector  (entire  structure) 
external  load  vector  for  an  element 

vectors  tangent  to  curvilinear  coordinate  lines 

(6  rj,  C) 

intensity  of  distributed  loads 

column  matrix  of  generalized  coordinates 

AAA 

displacement  vector  (=  ui  +  vj  +  wk) 

displacement  vector  of  node 

rotation  matrix  for  coordinate  transformation, 
Eq.  (3.40) 

area  of  curved  surface 
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U,  V,  w 


U  V  w  ■ 

t’  j’  ^ 


\U\ 

\Ua\ 


V. 

—  I 


His 


AAA 

V  V  V 
vV  v  V  V  3 


AAA 


^3* 


Vol 


w , 


wi 


y>  2 

x',y',z' 

xi > y,-*  2i 
a,  /3 
an  Pi 


y xy"1  yy2>  z 
^/A:,y,,  y  y'z"'>  l/;t:/z'’ 


shell  thickness  at  node  £ 

components  of  displacement  in  the  direction  of  *,  y  and  2 
axes,  respectively;  displacements  are  positive  in  the 
positive  direction  of  coordinate  axes 

components  of  displacement  at  node  i 

vector  of  nodal  parameters  for  entire  structure,  Eq.  (2.13) 
nodal  displacement  vector  for  an  element 
vector  of  parameters  at  node  i 

vector  of  parameters  at  node  i  for  shell  element, 

Eq.  (3.51) 

unit  vectors  in  directions  of  local  axes  x%  y '  and  2', 
respectively 

local  unit  vectors  at  node  i 

shell  thickness  vector  at  node  i 

volume  of  a  given  solid  domain 
work  done  by  external  load 

internal  work  of  strain  energy 

global  system  of  rectangular  coordinates 

local  system  of  rectangular  coordinates  for  shell  element 

coordinates  at  node  i 

rotations  of  nodal  normal  about  two  orthogonal  axes 
rotations  of  normal  at  node  i 

column  matrix  of  constant  coefficients,  Eq.  (3.21) 
shearing  strain  components 

shearing  strain  components  in  the  local  rectangular 
coordinates 
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8 


SUi 

U'] 


c 


C  ??)  f 

Vp 


[6] 


V 


[a] 


T  /  / ,  T  /  /,  T  /  / 
x y  ’  y  2  ’  zx 


[<£(«>  y>  s)] 
2 


L 


prefex  denoting  first  variation  of  a  function 
virtual  displacement  at  node  i  for  an  element 

local  strain  tensor,  Eq.  (3.42) 

local  strain  vector  for  a  shell  element  Eq.  (3.52) 

a  curvilinear  coordinate  in  the  thickness  direction  in  case 
of  a  shell  element 

curvilinear  coordinates  at  any  point  within  an  element 
curvilinear  coordinates  at  node  i 

position  constants  for  the  Gaussian  quadrature  point  i 

direction  casine  matrix  of  a  local  orthogonal  system  of 
axes 

Poisson’s  ratio 
stress  tensor,  Eq.  (3.43) 

shearing  stress  components  in  the  local  rectangular 
coordinates 

row  matrix  of  monomial  functions  of  cartesian  coordinate, 
Eq.  (3.1) 

summation  on  the  running  index  i 
integration  over  a  domain  D 
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CHAPTER  2 


THE  FINITE  ELEMENT  METHOD  OF  STRUCTURAL  ANALYSIS 
2.1  Background 

The  finite  element  methods  of  structural  mechanics  rely  heavily  on  numerical  compu¬ 
tation  and  their  advent  followed  the  availability  of  high-speed  digital  computers.  The  matrix 
formulation  of  structural  problems  was  formally  introduced  by  J.  H.  Argyris  in  a  series  of 
papers  published  in  1954—55.  About  the  same  period,  notable  progress  in  applying  finite  ele¬ 
ment  methods  to  the  analysis  of  aircraft  and  civil  engineering  structures  had  been  made  inde¬ 
pendently  by  a  number  of  investigators  including  M.  J.  Turner  (1956)  and  his  group  at  the 
Boeing  Company  and  B.  Langefors  (1958)  in  Sweden. 

The  basic  concept  of  the  finite  element  method  is  that  a  real  continuum  can  be  treated 
analytically  by  subdividing  it  into  a  finite  number  of  regions.  In  each  of  these  regions,  the 
behavior,  such  as  displacement  or  stress,  is  described  by  a  separate  field.  These  fields  are 
often  chosen  in  a  form  that  ensures  continuity  of  the  described  behavior  throughout  the  com¬ 
plete  continuum.  In  other  cases,  the  chosen  fields  do  not  ensure  continuity  but  nevertheless 
they  achieve  satisfactory  solutions.  These  later  cases  do  not  have  the  assurance  of  converg¬ 
ence*  possessed  by  the  fully  continuous  analytical  models  (Melosh,  1963  and  Irons  and  Draper, 
1965).  The  concept  of  the  finite  element  representation  owes  much  to  the  early  work  of 
A.  P.  Hrenikoff  (1941)  and  R.  Courant  (1943);  the  later  was  concerned  with  problems  governed 
by  broader  field  equations  than  just  structural  mechanics. 

Much  progress  had  been  made  in  the  various  aspects  of  finite  element  analysis.  Im¬ 
proved  results  were  realized  by  introducing  new  types  of  elements,  such  as  more  powerful 
refined  elements  (Argyris,  1965;  Felippa,  1966;  Mehrain,  1967;  Kohnke  and  Schnobrich, 

1969;  and  Chu  and  Schnobrick,  1970)  or  efficient  superelements  (Araldsen  and  Egeland, 

1971).  Successful  developments  were  also  cited  for  various  forms  of  structural  behavior 
representation  as  in  dynamics,  plasticity,  and  large  deflection  (Argyris,  1965;  Przemieniecki 
et  al.,  1971;  and  Zienkiewicz,  1971). 

The  formulation  of  the  finite  element  method  can  be  traced  to  energy  procedures,  prin- 
cipally  the  minimum  potential  energy  (MPE)  method  and  the  minimum  complementary  energy 
(MCE)  method.  The  MPE  method  is  associated  with  assumed  displacement  parameters  as  un¬ 
knowns  and  is  usually  termed  the  “displacement”  or  “stiffness”  method.  On  the  other  hand, 
the  MCE  method  deals  with  stress  parameters  and  is  termed  the  “flexibility”  approach.  Sev¬ 
eral  authors,  e.g.,  Fraeijs  deVeubeke  (1964),  have  been  concerned  with  a  parallel  use  of  dis¬ 
placement  and  equilibrium  models  to  obtain  lower  and  upper  bounds  to  the  exact  solution. 

Still  others  (Fraeijs  deVeubeke,  1964  and  Herman,  1967)  used  a  mixed  model  and  considered 
both  displacements  and  stresses  as  primary  variables.  The  ease  with  which  a  continuous 


*For  convergence  requirements,  see  Section  2.3.1. 
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displacement  pattern  can  be  prescribed  (compared  to  the  alternative  approach  of  forming  an 
equilibrating  internal  force  field)  has  aided  the  widespread  use  and  development  of  the  finite 
element  displacement  approach.  The  displacement  model  and  the  stiffness  analysis  are  em¬ 
ployed  in  the  present  study. 

2.2  Finite  Element  Displacement  Approach 

The  displacement  formulation  involves  derivation  of  the  stiffness  matrix  of  each  indi¬ 
vidual  element.  The  stiffness  matrix  of  the  entire  assembled  structure  is  then  obtained  by 
the  direct  stiffness  method.  This  matrix,  along  with  the  prescribed  displacement  boundary 
conditions  and  loads,  is  used  for  the  solution  of  displacements  and  stresses. 

2.2.1  Element  Analysis 

The  basic  steps  for  derivation  of  the  element  stiffness  matrix  are: 

a.  Express  the  internal  displacements  \U\  of  the  element  in  terms  of  displacement  func¬ 
tions  96  (rr,  y,  3),  and  generalized  coordinates  \q\ 

Wi  =  [<£]{?!  (2.1) 

where  tu  ( x ,  y,  3)' 

1C}  =  v  ( x ,  y,  z)  ■ 
w{x,  y,  2), 

is  a  displacement  vector  consisting  of  displacement  components  u,  v,  w,  referenced  to  the 
rectangular  cartesian  coordinate  system  (x,  y,  z). 

b.  Express  the  nodal  displacements  \Ua\  in  terms  of  generalized  coordinates  Igd: 


\Ua\  =  [A]\q\  (2.2) 

Here  \Va\  =  {t/(.(a r,  y,  z)\  and 

Ui(x,  y,  2)  =  V  (*.,  yv  2.)  i  =  1,  2,  . .  ,  n 

are  the  displacements  at  node  and  n  is  the  number  of  nodes.  Coefficients  of  matrix  [/l] 
are  functions  of  nodal  coordinates  (a:.,  y .,  3^).  Conversely, 

ijuwr1  n  (2.3) 
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c.  Evaluate  strains  {«!  from  the  assumed  displacement.  Use  the  strain-displacement 
relationship 


U\  =  £{£/}=  lG(x,  y,  a)]  \q\ 


=  [G]  [A-1]  \Va\ 


or 

M=[S]{t/“l  (2.4) 

where  £  is  a  differential  operator  and  [<?]  =  £  [<£]. 

d.  Compute  stresses  jai  using  the  elasticity  matrix  [D]  established  from  the  properties  of 
the  material 


M  =  [D]  U! 


(2.5) 


or 


M  =  [D]  [B(x,  y,  a)]  \Ua\  (2.6) 

e.  The  condition  for  equilibrium  is  obtained  by  applying  the  principle  of  minimum  potential 
energy 


d(wi  +  w  e )  =  0 

During  a  virtual  displacement  Ua  ->  SUa,  the  internal  work  done  is 

8Wi  =  8fU \T  iaU(Vol) 

=  [D]  [ B ]  c?(Vol)  {Ua\ 


(2.7) 


or 

8v>.  =  \8Ua}T  [Ke]  \Ua\  (2.8) 

The  external  work  done  by  a  set  of  nodal  forces  \Pa]  corresponding  to  the  nodal  displacements 
{Ua\ is 

8we  =  -  \8Ua]T  \Pa\  (2.9) 


Substitute  into  Eq.  (2.7).  Since  the  virtual  displacements  8Ua  are  arbitrary,  we  have 


where 


[/g  =  f[B]r[D][B]  d(V ol) 

=  u-l]T $mT mm  <*(v 0i)  ur1]  (2.11) 

expresses  the  nodal  force-displacement  relation  and  is  the  desired  element  stiffness  matrix. 

f.  Establish  the  load  generalization.  Generally  loads  are  distributed.  Concentrated 
loads  at  nodes  represent  special  cases.  In  the  finite  element  context  where  all  the  forces 
can  be  transmitted  only  through  the  prescribed  network  of  nodes,  we  need  to  compute  the 
equivalent  concentrated  nodal  loads  f Pa\  for  the  actual  distributed  loading  p  (x,  y,  g). 
Equivalence  is  based  on  the  work  done  during  a  virtual  displacement  consistent  with  the 
assumed  displacement  field  U  (x,  y,  z).  Since 

8we  =  -SfW\T  {pi  d  (Vol) 

Eq.  (2.9)  gives  the  equivalent  loading  vector 

iPj  i  =  LT1]7/^]7,  {pi  d (Vol) 
or 

!P;i=/[N]rlpU(Vol)  (2.12) 

The  vector  S P^  \  is  often  called  the  generalized  load  of  element  “a”  and  N  is  known  as  shape 
function.  Thus  a  normal  load  can  produce  not  only  parallel  nodal  forces  but  also  couples,  or 
their  equivalent,  and  these  will  depend  on  the  displacement  assumption  U  (x,  y,  z)  used. 

2.2.2  Structural  Analysis  (by  Direct  Stiffness  Method) 

The  real  elastic  structure  is  now  represented  by  a  finite  number  of  small,  discrete 
elements.  Once  their  approximation  behaviors,  identified  by  their  individual  stiffness  matri¬ 
ces  [/f“],  have  been  established,  the  stiffness  matrix  [K]  for  the  complete  structure  is  ob¬ 
tained  by  the  proper  summation  of  each  element  stiffness  matrix  in  the  structure.  This  is 
done  conceptually  by  joining  successive  elements  at  their  adjoining  nodes  and  requiring  that 
the  conditions  of  compatibility  of  displacements  and  equilibrium  of  forces  have  been  satisfied 
at  every  node  throughout  the  structure.  A  set  of  simultaneous  equations  is  generated  in  terms 
of  displacement  parameters  j U\.  For  compatibility,  express  in  matrix  form: 

\Ua\  =  [La]  \U\  (2.13) 
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where  \Ua  \  is  the  column  vector  of  element  nodal  displacements, 

\U\  is  the  column  vector  of  nodal  displacements  for  the  entire  structure, 

[La]  is  the  localizing  matrix  for  element  “a”  relating  the  displacements  of  the  element 
to  those  of  the  structure  and  is  with  respect  to  a  set  of  nodal  numbering  system, 
and  for  equilibrium, 


[K]  It/j  =  |P  |  (2.14) 

M  T 

[K]  =  2  [La]  [/7“]  [La] 

0=1 

is  the  structure  stiffness  matrix, 
is  the  element  stiffness  matrix, 
is  the  total  number  of  element  in  the  structure,  and 
is  the  structure  loading  vector  and  is  equal  to 

M  T 
{P\=  1  [La ] 

a=  1 

assembly  process  is  also  known  as  the  direct  stiffness  method.  Matrix  notation 
is  convenient,  general,  and  applicable  to  a  wide  range  of  structural  problems.  In  practice, 
further  processing  takes  into  consideration  the  fact  that  the  structure  stiffness  matrix  is 
banded,  symmetric,  and  sparsely  populated.  This  allows  a  significant  advance  (Tezcan, 

1966;  Irons,  1970;  and  Jennings  and  Tuff,  1971)  in  computational  efficiency  (the  number  of 
arithmetic  operations  and  data  storage  can  often  be  drastically  reduced). 

2.3  Characteristics  of  Finite  Element  Analysis 
2.3.1  Convergence  Criteria 

The  trace  of  element  characteristics  follows  a  prescribed  path  once  the  shape  function 
N  has  been  defined.  It  describes  the  displacement  field  within  the  element  in  terms  of  the 
nodal  parameters.  The  shape  function  matrix  [/V]  can  be  determined  from  Eq.  (2.1),  (2.3),  and 
(2.12): 


Note  that 


where  177] 

[i Ka ] 

M 

\P\ 


The 


and 


!£/!  =  [/V]  \Ui | 

(2.15) 

[/V]  =  [0]  [A]-1 

(2.16) 
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The  shape  functions  are  usually  taken  as  polynomial  expressions  in  terms  of  either  the  global 
coordinates  (%,  y,  z)  or  the  local  coordinate  system  (f,  77,  £). 

The  reliability  of  a  solution  by  the  finite  element  method  is  indicated  by  the  fact  that 
the  approximate  numerical  results  come  increasingly  closer  to  the  correct  value  as  the  finite 
element  mesh  is  repeatedly  subdivided  into  finer  and  finer  meshes  (Tong  and  Pian,  1967  and 
deAranetes,  1968). 

The  requirements  of  convergence  fall  into  two  categories:  (a)  completeness  of  the  dis¬ 
placement  field  V  and  (b)  interelement  continuity.  First,  the  completeness  requirement  en¬ 
sures  that  the  energy  represented  by  the  functional  (such  as  the  integral  of  potential  energy) 
includes  a  constant  energy  state  for  each  element.  If  this  is  satisfied,  the  true  energy  state 
of  the  entire  structure  can  be  represented,  in  the  limit,  as  the  mesh  layout  is  refined.  Mathe¬ 
matically,  this  requires  that  all  uniform  states  of  the  displacement  variables  V  must  be  in¬ 
cluded  in  each  element.  This  results  in  the  equivalent  requirement  that  “constant  strain” 
and  “rigid  body”  states  must  be  included. 

The  second  requirement  is  that  of  continuity  between  adjacent  elements.  In  order  for 
the  functional  at  the  element  interface  to  remain  finite,  it  has  been  considered  necessary  to 
provide  continuity  of  displacement  variable  U  and  its  derivatives  to  one  order  lower  than  the 
order  of  the  highest  derivative  of  that  variable  in  the  energy  integral.  The  success  of  certain 
nonconforming  elements,  however,  has  led  to  a  reevaluation  of  this  requirement.  Bazeley  et  al. 
(1965)  suggested  a  less  stringent  condition,  namely,  that  reduced  continuity  must  be  maintained 
for  the  state  of  constant  energy  in  the  region  concerned. 

2.3.2  Elements  of  Arbitrary  Shapes 

Linear  elasticity  problems  can  be  readily  solved  by  a  finite  element  technique  that  em¬ 
ploys  simple  triangular  or  tetrahedral  elements.  The  displacement  fields  for  those  elements 
have  often  taken  the  form  of  polynomials  of  variables  in  cartesian  coordinates.  That  is 

\U\  =  [<f>{x,  y,  z)]  \q\ 

or 

u  (x,y,z)  =  gn  +  q12x  +  q13y  +  q14z  +  .  .  .  \ 

v  (x,y,z)  ^  q21  +  q22x  +  q23y  +  q24z  +  .  .  .  \  (2.17) 

w(x,y,z)  =  q31  +  q32x  +  q33y  +  q34z  +  .  .  .  / 

In  the  case  of  the  constant  strain  element,  the  criteria  of  convergence  can  be  satisfied  when 

the  assumed  displacement  functions  include  the  complete  first  order  polynomial  [1,  *,  y,  s] 
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in  addition  to  the  basic  requirement  of  compatibility*  of  nodal  parameters.  Additional  nodal 
variables  can  be  introduced  to  give  better  element  performance  with  improved  representation 
of  the  actual  deformation.  Selective  higher  order  terms  are  to  be  added  to  ensure  the  continu¬ 
ity  of  displacement  V  across  the  interelement  boundary. 

The  popularity  of  polynomial  displacement  expansions  </>(#,  y,  z)  lies  in  the  fact  that 
matrix  operation  is  straightforward  and  stiffness  coefficients  can  be  readily  obtained  in  ex¬ 
plicit  form  by  a  standard  procedure  (see  Section  2.2).  It  becomes  obvious  in  application,  how¬ 
ever,  that  the  property  of  a  displacement  variable  U  (x,  y,  z)  should  have  no  preferred  directions. 
This  requires  that  the  complete  polynomial  expansion  be  used  in  the  displacement  assumption. 
To  implement  this,  a  simple  procedure  can  be  found  to  determine  shape  functions  [ N ]  and  dis¬ 
placement  patterns  \U  i;  the  tacit  assumption  is  that  matrix  [A]  (Eq.  (2.14))  is  invertible  and 
that  the  shape  function  is  correct  simply  because  there  is  a  “match”  between  the  monomials 
present  in  the  displacement  expansion  and  the  corresponding  nodal  variables  (Dunne,  1968). 
However,  this  approach  is  not  generally  valid. 

Complete  polynomials  for  a  compatible  displacement  field  are  well  suited  for  the  tri¬ 
angular  and  tetrahedral  family  of  elements,  but  the  use  of  complete  polynomial  displacement 
expansion  may  lead  to  algebraic  difficulties  in  the  case  of  arbitrary  quadrilaterals,  hexahedra, 
and  plate  bending  elements.  On  the  other  hand,  it  is  possible  to  formulate  a  compatible  ele¬ 
ment  that  possesses  rotational  invariance  via  direct  formulation  by  employing  a  coordinate 
transformation  (mapping)  or  a  natural  coordinate  system  (Irons,  1966  and  Ergatoudis  et  al., 

1968).  Additional  kinematic  capabilities  can  be  incorporated  by  introducing  intermediate  nodal 
parameters  (Zienkiewicz  et  al.,  1971).  This  is  another  desirable  feature  of  isoparometric 
elements.  A  more  detailed  description  will  be  given  in  the  following  chapters. 


*This  ensures  the  minimum  continuity  requirements  of  displacement  at  element  interface. 
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CHAPTER  3 


ELASTIC  ANALYSIS  IN  THREE-DIMENSIONAL  SPACE 

3.1  Introduction  to  Solid  Elements 

For  many  years,  considerable  effort  has  been  devoted  to  solving  problems  in  the  realm 
of  three-dimensional  elasticity.  The  classical,  analytical  approach  for  solving  a  set  of 
governing  differential  equations  derived  from  a  three-dimensional  theory  is  available  only  for 
bodies  with  simple  geometric  forms  and  for  restricted  boundary  conditions  and  limited  loading 
cases.  Different  numerical  procedures  (finite  difference  methods  etc.)  have  been  applied  to 
solve  these  differential  equations;  their  success  is  generally  still  limited  to  special  geometric 
shapes  and,  at  most,  is  of  occasional  academic  interest.  Now  that  the  finite  element  method 
has  been  eminently  successful  in  dealing  with  certain  complex  problems  such  as  plane  stress 
and  plate  bending,  it  has  perhaps  even  greater  potential  for  the  solution  of  three-dimensional 
problems. 

The  approach  envisioned  here  is  based  on  a  full  three-dimensional  analysis  and  is  com¬ 
pletely  general  in  nature.  It  will  be  shown  that  the  solid  elements  used  in  the  present  study 
are  capable  of  correctly  representing  the  behavior  of  a  beam,  plate,  shell,  or  any  of  the  varied 
aspects  of  structural  components. 

Because  the  tetrahedron-shaped  element  (Fig.  3.1)  has  simplicity  and  flexibility,  it 
was  a  natural  choice  for  the  earlier  development  of  solid  elements  (Gallagher  et  al.,  1962). 

The  element  shape  is  defined  by  four  arbitrary,  noncoplanar  points  in  space.  Any  topography 
may  therefore  be  represented  with  sufficient  accuracy  by  some  assembly  of  these  tetrahedron 
elements.  The  drawback  to  this  element  shape  is  the  large  number  of  element  inputs  required 
to  describe  a  complex  surface.  The  mesh  is  difficult  to  visualize,  and  the  amount  of  data  to 
be  processed  as  well  as  the  number  and  the  bandwidth  of  the  system  of  equations  generated 
tend  to  exceed  the  storage  capacity  of  the  average  size  computer  and  call  for  excessive  com¬ 
puter  times.  To  reduce  these  constraints,  isoparametric  elements,  including  those  with  a 
curved  face,  have  been  introduced.  These  elements  represent  a  great  improvement  over  the 
tetrahedron  because  they  enable  bodies  with  curved  boundaries  to  be  treated  with  a  limited 
number  of  elements. 

3.2  The  Basic  Solid  Element 

A  comparison  (Clough,  1969)  of  the  performance  of  solid  finite  elements  has  shown  that 
isoparametric  hexahedron  elements  are  distinctly  superior  to  any  tetrahedron  assemblages,  both 
in  terms  of  the  properties  of  the  individual  elements  and  their  application  to  analyze  real  struc¬ 
tural  systems.  Isoparametric  elements  have  the  additional  advantage  of  isotropy.  It  is  evident, 
then,  that  a  general-purpose,  three-dimensional  finite  element  analysis  program  should  be  con¬ 
structed  around  the  isoparametric  hexahedron  element  family. 
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3.2.1  The  Isoparametric  Displocement  Field 


Consider  the  general  family  of  elements  with  six  faces,  the  hexahedron  family.  The 
use  of  an  appropriate  natural  coordinate  system  greatly  simplifies  the  formulation  of  the  ele¬ 
ment  stiffness  matrix  for  some  very  complex  members  of  this  family.  These  coordinate  lines 
are  generally  curved  in  space  and  follow  only  the  interface  topology  of  an  individual  element. 
As  shown  in  Fig.  3.2  (and  later  in  Fig.  3.3)  for  instance,  f,  y,  C  are  natural  coordinates;  each 
coordinate  axis  is  associated  with  a  pair  of  opposing  faces  which  are  given  the  coordinate 
values  of  +  1.  In  their  local  reference  frame,  the  elements  take  on  the  image  of  a  2  by  2  by  2 
cube,  whereas  in  the  real  cartesian  coordinate  system  ( x ,  y,  z),  they  can  be  any  arbitrarily 
warped,  six-sided  solids. 

We  begin  with  the  linear  hexahedron  that  has  the  eight  corner  nodes  shown  in  Fig.  3.2 
and  introduce  a  polynomial  expansion  of  the  displacement  component  in  rectangular  cartesian 
coordinates 

u  =  u  (x,  y,  z) 

=  ?n  +  ?i2r+  ?i3  y  +  ?i42  +  ?i5  xv +  Qif>yz  +  ?i7 zx  +  vi&xvz 
or 

u  =  [<f>i(x,y,z)]  i=l,2,  ....,8  (3.1) 

It  will  be  demonstrated  (Section  3.2.4)  that  the  criteria  of  convergence  are  satisfied.  Without 
losing  generality,  we  can  consider  an  element  in  the  form  of  a  rectangular  parallelopiped  with 
the  origin  located  at  centroid  and  the  coordinate  planes  as  planes  of  symmetry.  By  evaluating 
nodal  displacement  parameters  at  the  appropriate  coordinates,  i.e.,  ui  =  M  (on,  y(.,'a  .),  etc.,  we 
can  assemble  the  nodal  displacement  vector,  for  instance,  luti  =  [X-]  etc.  From  Eq. 

(2.15),  we  obtain  by  identity 


(3.2) 

,  8  (3.3) 

where  xi  =  ±  1,  y.  =  ±  1,  and  s  .-=  + 1  are  coordinate  values  at  node  i. 

Here,  N.  ( x ,  y,  z)  is  known  as  the  shape  function,  or  interpolation  function,  such  that  it  takes 
on  unit  value  at  the  indexing  node  ( i )  and  is  zero  at  all  other  nodes.  These  interpolation  func¬ 
tions  are  readily  developed  simply  by  writing  them  as  the  product  of  the  equations  of  the  lines 
or  surfaces  through  all  the  remaining  nodes.  This  can  frequently  be  done  by  direct  inspection. 


u  =  Nlu1  +  N2u2  + . +  /V8  us 

8 

=  2  N.u. 

»=  l  *  ‘ 


and 


Ni  =  ~  (!  +  xix)  (i  +  yty)  (!  +  st  s)  *  =  i, 
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Figure  3.2  An  eight-node  generalized  hexahedron 


w 


Figure  3.3a  Quadratic  element  (20  nodes) 


LOCAL  COORDINATES 


Figure  3.3b  Cubic  element  (32  nodes) 

Figure  3.3  Refined  curved  hexahedron 


as 


Proceeding  in  the  same  fashion,  we  can  define  other  displacement  components  v  and  w 


v  =  v(x,  y,  z) 

8  "I 

=  2  AC  v.  I 

: 

w  =  £  N .  w  ■  i  =  1,2, . ,8  /  (3.4) 

i=  l  1  1 


where  the  shape  function  AC  takes  on  the  same  form,  Eq.  (3.3). 

As  mentioned  earlier,  an  element  of  prismatic  shape  has  rather  limited  appeal  in  deal¬ 
ing  with  practical  problems.  This  geometric  constraint  can  be  relaxed  by  a  suitable  coordinate 
transformation  commonly  known  as  mapping.  Consider,  for  instance,  an  element  of  Fig.  3.2a 
being  distorted  geometrically  to  a  shape  shown  in  Fig.  3.2b.  In  other  words,  map  the  element 
into  x,  y,  2-coordinate  space  such  that  a  typical  node  i  moves  to  a  position  (x.,  y z^.  A 
relation  of  the  forms 

x  =  1  N .  x. 

i  1  1 

y  =  2  N.  y- 

2  =  1  AC  si  i  =  1,  2,  .  .  .  ,  n 


(where  n  -  NNPE,  the  number  of  nodes  per  element)  gives  »  =  #(£,  rj,  £),  .  .  .  and  can 
be  used  to  define  the  mapping.  Here,  the  shape  functions 

AC  =  A',  (£  y,  C) 

=  -  +  +  +  (3.6) 


are  written  in  terms  of  the  local  dimensionless  coordinates  £,  y,  and  £.  Following  the  nodal 
ordering  number  as  given  in  Fig.  3.2a,  we  have  the  local  nodal  coordinates 

(£,  Vrti)  =  (-!,-  !,-!),(  1,-1,  -4), 

(-1,  !,-!),(  1,  1,-1), 

(-1,-1,  1),  (  4,-4,  4), 

(-1,  1,  1),  (  1,  4,  1), 
for  i  =  1,  2,  .  .  .  ,  8 
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The  displacements  u,  v,  w  can  be  interpolated  in  the  same  manner;  hence  u  =  u  (£,  -q,  £), 
etc.  Again,  we  write 


u  =  2  N .  u.  1 

i  1  1  I 

v  =  2  N .  v.  t 

l  l  I 

w  =  2Ni  w.  i  =  1,  .  .  .  .  ,  n  / 
i 


(3.7) 


The  shape  functions  7/^,  used  here  are  identical  to  those  employed  in  the  coordinate 

definition  and  in  such  cases  the  formulation  has  been  termed  “isoparametric”  (Zienkiewicz 
et  al.,  1969).  With  the  position  and  the  displacement  of  a  material  point  in  an  element  space 
defined,  a  family  of  isoparametric  elements  can  be  formed  routinely. 


3.2.2  Numerical  Calculation  of  Stiffness  Matrix 

Consider  an  element  stiffness  matrix,  Eq.  (2.11), 


[Kel  =  fff  W]1[D][B]  dxdydz 

volume  of 
element 


To  integrate  it,  we  need  to  evaluate  the  integrand  [£]r[£>][B]  =  G(x,  y,  z)  in  terms  of  independ¬ 
ent  space  variables  (x,  y,  z).  The  algebraic  expression  of  G  involves  the  strain  matrix  [ B ] 
which  is  derived  from  the  definition  of  strain  and  consists  of  appropriate  first  derivatives  of 
the  displacements.  For  a  solid  element,  we  have 


X 


u  =  < 


xy 


y z 


Y  Z  X 


dN. 

1 

dx 

0 

0 


dy 

0 

dNi 

dz 


dN. 

1 

dy 

0 


dN.  dN- 


dz 

0 


dN. 

I 

dz 


dx 

dN.  dN. 


dy 

dN. 

I 

dx 


w . 


*'=M, 


(3.8) 


NNPE 


[B]  {Ua\ 

6  x  24  24  x  1 
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for  NNPE  =  8 


It  is  clear  then  the  expression  G  is  algebraicly  complex  and  its  integration  presents  a  formid¬ 
able  task  since  the  shape  functions  N  i  =  N 77,  £)  are  established  in  local  rather  than  global 
coordinates.  In  order  to  proceed  further,  it  is  desirable  to  introduce  an  auxiliary  expression 
g(£,  tj ,  C)  in  terms  of  local  coordinates  (f,  77,  £). 

The  Jacobian  of  the  transformation  is  defined  as  J  =  d(x,  y,  s)/d{£,  77,  £)  and  the 
matrix  is 


dN. 

where  - 

H 

d/v. 

d  rj 

dNi 

H 


dx 

dy  da 

d|  di 

d  a; 

dy  d  2 

.  d  77 

d  77  d  77 

dx 

dy  d  2 

Jc 

I? 

1  VP 

1  ^ 

Hn1 

dtt2 

dNg 

d£ 

d£ 

d{ 

dN  j 

dN  2 

00 

<"C> 

d?7 

d  V 

dr] 

d^i 

dN2 

dN& 

_d<T 

dC 

d£_ 

it 

=  (1  +  ViV)  (!  +  CiC) 

(1  +  £■£) 

=  — —  0i,)(W,O 

=  (i  +  £,£)  (i  +77777)—  ,  1  =  1,  2, 


,  8 


^8 


(3.9) 


Following  the  standard  rule  of  partial  differentiation,  for  example, 
dN;  dN:  k  d/V..  d/V. 


da;  d  £  da;  d  77  da; 


d<T 


d£  d: 


,  etc., 
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we  have 


dNt 

dN  2 

00 

d  x 

dx 

d  x 

<9A1 

dN2 

00 

dy 

dy 

dy 

dNx 

^  2 

00 

d  z 

d  2 

d  2 

dNx 

dN  2 

d  N  q 

d£ 

<9/V1 

dN  2 

00 

d-q 

dy 

"'dr, 

dNx 

dN  2 

00 

K 

d  C 

■  "  dC 

(3.10) 


Now  the  values  of  dN  t/d^,  dN  x/ dr,,  dN  t/dCi  and  matrix  [=7]  can  be  calculated  point  by 
point  in  an  element  subregion.  The  matrix  [B\  can  be  readily  assembled.  The  components  of 
strains  everywhere  within  the  elastic  element  are  now  defined  in  terms  of  nodal  deflections 
\Ua\  as  parameters. 

The  integration  in  terms  of  local  coordinate  variables  (f,  rj,  C)  can  be  executed  in  a 
simple  format,  that  is, 


[K ]  =  fff  G(x,  y,  z)  dx  dy  dz 
e  Vol 

=  [B]\J\d£dvdC  (3.11) 

fftf,  V,  O  =  [B]T[D]  [B]  \J({,  r,,  C) I  (3.12) 


For  all  practical  purposes  the  integration  performed  numerically  (Hammer,  1959)  as  a 
summation  of  quadratures  is  a  simple  matter  on  a  digital  computer.  Here, 


».]“ 


NPZ  NPY  NPX 

2  2  2 

iz— 1  iy=l  ix= 1 


’ll ..  (i,)H 


H  H 

ix  n  iy  n  iz 


(3.13) 


If  the  Legendre-Gauss  quadrature  formula  is  employed,  tjix,  ? ^  ,  (iz  correspond  to  the  abscissas 
for  the  zeros  of  the  Legendre  polynomial  of  degree  NPX,  .  .  .  etc.  NPX,  NPY,  NPZ  are  the 
number  of  integration  points  to  be  used  in  each  linear  quadrature  space.  The  values  of 
and  its  weight  factors  Hix  can  be  found  in  standard  quadrature  tables  (Stroud  and  Secrest, 

1966). 

It  is  worth  noting  that  the  stiffness  matrix  [ K  ]  is  symmetric  and  that  the  strain  matrix 
[6]  and  the  elasticity  matrix  [Zl]  include  many  zero  terms,  or  null  submatrices.  A  substantial 
reduction  in  arithematic  operation  is  possible  by  carrying  the  algebraic  processing  a  little 
further.  First,  the  strain  matrix  [B] 
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[. B ] 

6  x  24 


(3.14) 


^21 

_f„_i 

R  1 
°22  1 

i 

^23 

~B\ 

0 

0~ 

0 

62 

0 

0 

0 

63 

62 

61 

0 

0 

63 

62 

_63 

0 

51_ 

can  be  expressed  in  terms  of  submatrices: 


61  (8)  = 

dNt 

dN2 

dN& 

dx 

’  dx  ’ 

dx 

62  (8)  = 

dNx 

dN2 

00 

<t> 

dy 

’  dy  ’ 

dy 

63  (8)  = 

~dN1 

dN2 

dN~ 

dz 

'  dz  ' 

’  dz 

(3.15) 


Second,  the  elasticity  matrix  [6]  which  relates  the  stresses  and  strains,  Eq.  (2.5),  can  include 
any  anisotropic  properties  and  can  be  prescribed  as  a  function  of  spatial  distribution,  e.g., 
y,  z).  The  [6]  matrix  is  symmetric  and  can  be  written  as: 


where 


6 

0 

[6(6,  6)]  = 

m 

o 

6 

LiJ 

■*» 

^12  ^13 

0J- 

SYM 

^22 

^23 

^33 

^44 

0 

0 

0.1- 

SYM 

DSS 

0 

°66_ 

(3.16) 


[0]  is  a  3  by  3  null  matrix.  For  this  study,  we  shall  consider  an  isotropic  material. 


that  case 
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CC 
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to.] -  OK 
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[6s]  =  6tf 
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SYM 
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SYM 
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where  DK  - 


E 


(1  +  v)  (1-  2v) 

AA  =  l-i/ 

BB  =  v 

CC  =  (1  -  2 1/)/2 

£  and  i/  are  the  usual  elastic  constants,  Young’s  modulus  and  Poisson’s  ratio,  respectively. 

At  this  stage,  the  calculation  of  element  stiffness  matrix  can  be  broken  into  parts. 
Terms  that  contain  null  factors  can  be  screened,  large  number  of  intermediate  calculations 
can  be  eliminated,  and  economy  in  computer  time  can  be  realized.  For  instance: 

[Ke  (i, ;)]  =  /I  [B  (tn ,i)]T  [D  (m,  n)]  [B  ( n , ;)]  \J\  d £  d n  d C  (3.17) 

24x24  24x  6  6x  6  6x24 


HU 


or 

i —  — 

K  ( r ,  s)  K  (r,  s  +  8)  K  ( r ,  s  +  16) 

£(r  +  8,  s  +  8)  £(r+8,  s  +  16)  (3.18) 

_SYMM  ^(r+16,  s  +  162 

r,  s  =  1,  2, . ,8 

A  typical  submatrix,  for  example,  is 

=  fiu  +  e27!  621]|j|  *  £<*,<*{  (3.19) 

This  can  be  reduced  to 

Dn  B1(s )  +  ffi2(0  ^44  S2(«)  +  ^3(r)  Z?66  63(s)j]  1=71  d{dr,dC  (3.20) 

Other  submatrices  are 


/'/'/MSiW  D12B2(s)  +  B2(r)  D 44  61(S)]|=/|  rff  tfr? 

i\  L\  L\  W1 W  °  13  fi3  +  63  <f>  ^66  fil  («)1  V  dC 

f_\[B2(r)  D22B2(s)  +  B\(r)  D44Bl(s)  +  B3(r)  D55  63(a)]  |=/|  <7  77 

/Xi  J1  J1  [B2(r)  D23  B3(s)  +  B3(r)  Z?55  B2(s)]  iJjd^d^dC 
f  f  /  [58(f)  Z?33  63  (s)  +  62 (r)  6g5  62(s)  +  61  (r)  Z?66  61(a)]  |7|  d^drjdC 
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For  submatrices  on  the  diagonal,  stiffness  coefficients  must  be  computed  only  for 
those  terms  where  s>r,  since  advantage  is  taken  of  the  symmetry  of  the  stiffness  matrix. 
Finally  the  nodal  parameters  can  be  regrouped,  i.e., 

Wa\  =  \U.\ 


W  ={  •/ 


i=  1,2, . ,8 


to  expedite  the  assembly  into  the  entire  structural  system  matrix  [K~\  of  Eq.  (2.14). 

3.2,3  Higher  Order  Curved  Elements 

In  earlier  developments,  the  boundary  of  a  finite  element  was  usually  thought  of  as 
composed  of  straight  edges  and  its  geometric  form  completely  defined  by  a  series  of  corner 
nodes.  Where  an  edge  is  curved,  appropriate  intermediate  nodes  must  be  specified  along  that 
edge.  One  intermediate  point  will  suffice  to  define  the  shape  of  an  edge  with  a  simple  or 
constant  curvature.  Two  or  more  intermediate  points  will  be  required  to  specify  the  geometry 
of  an  edge  with  a  multiple  or  reversed  curvature. 

3.2.3. 1  Quadratic  Curved  Element.  In  most  cases,  the  boundary  of  a  complex  structural 
shape  can  be  closely  represented  by  a  sequence  of  quadratic  curves  or  quadratic  surfaces. 
Therefore  many  complicated  problems  can  be  realistically  defined  and  solved  with  only  a 
limited  number  of  simple,  discrete,  curved  elements.  These  curved  elements  can  be  readily 
formulated  by  using  isoparametric  concepts  (Section  3.2.1).  For  example,  by  placing  one 
midside  node  along  each  edge  of  an  8-node  hexahedron,  we  obtain  a  20-node  element.  This 
20-node  hexahedron  is  coming  into  widespread  use  and  is  an  important  tool  in  three- 
dimensional  analysis  (Fig.  3.3a). 

For  a  20-node  hexahedron  element  with  edges  capable  of  displaying  simple  curvatures, 
we  can  write  a  displacement  expansion  that  contains  12  terms  in  addition  to  those  required 
for  the  8-node  hexahedron  (Eq.  (3.1))  to  ensure  compatibility. 

Hence, 

w20  =  [<£;(*>  SS  K-!  *  =  1,2, . ,20 

=  «j  +  a2%  +  a3y  +  a4z  +  a^xy  +  a6yz  +  a^zx  +  a&xyz 
+  agx2  +  a10y2  +  aua2 

+  al2x2y  +  a13x2z  +  a14y2x  +  a15y2z  +  a16z2y  +  alyz2x 


+  +  al9xy2z  +  a^xyz2 


(3.21) 
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It  is  seen  that  the  displacement  is  capable  of  a  quadratic  variation  along  any  space  variable, 
and  consequently  the  20-node  element  is  often  called  a  quadratic  element. 

Following  the  example  of  Section  3.2.1,  we  obtain  immediately  the  shape  functions  for 
the  20-node  hexahedron  element: 

For  corner  nodes  i  =  1  to  8,  =  ±  1,  n.  =  ±  1,  £,•  =  ±  1 

Ni=  ~Q-  +  €ie)(l+'liri)V-+CiO(€i€+Viri  + CiC-2)  (3-22) 

For  midside  nodes  i  -  9  to  12,  f  ■  =  +  1,  r)i  =  +  1,  £.  =  0 

/v.=  j(i  +  <f.O(i  +  ^)(W2) 

and  midside  nodes  i  =  13  to  16,  ■  =  +  1,  =  0,  £.  =  ±  1 

w.=  Y(i  +  ^.0(i->/2)  (i  +  ^-O 

and  midside  nodes  i  =  17  to  20,  =  0,  r\i  =  +  1,  £{.  =  +  1 

^•  =  j(i-  ^Mi  +  if ^)(i+ 

The  displacements  and  coordinates  definition  for  the  quadratic  element  follow  directly 
Eq.  (3.7)  and  (3.5)  of  the  previous  case.  The  summation  will  extend  to  all  20  nodes,  i.e. 

Iu  =  *  Nt  u. 

v  =  2  N .  v  . 
i  1  1 

w  =  2  N i  w i  i=  1,  2, . ,20 

The  calculation  of  element  stiffness  matrix  follows  the  procedure  outlined  in  Section 
3.2.2.  For  the  quadratic  element,  matrix  [/f6]  is  of  an  order  of  60  x  60;  a  Fortran  listing  of  a 
subroutine  for  computing  stiffness  coefficients  for  this  element  is  included  in  the  appendix. 


(3.23) 
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3. 2. 3. 2  Cubic  Curved  Element.  By  placing  two  intermediate  nodes  at  the  one-third  and  two- 
third  points  along  each  edge  of  an  8-node  hexahedron,  we  then  obtain  a  32-node  cubic  element, 
Fig.  3.3.  To  ensure  compatibility,  a  displacement  expansion  involving  32  terms  can  be 
written: 

« 32  (®>  s)  =  U2o(X,  y ’  2 )  +  a21x3  +  a22^3  +  a23s3 

+  a24x3y  +  a25x3z  +  a26y3z  +  a2yy3x  +  a2gz3x  +  a2gz3y 

+  a3Qx3yz  +  a31xy3z  +  a32xyz3  (3.24) 

where  the  first  20  terms  are  identical  to  u2Q  in  Eq.  (3.21). 

The  32  shape  functions  for  the  cubic  element  can  be  obtained  by  standard  procedure  or 
by  inspection  (Zienkiewicz,  1971)  and  the  96  x  96  element  stiffness  matrix  can  be  computed 
in  a  manner  similar  to  those  already  described. 

3.2.4  Practical  Considerations 

It  is  appropriate  to  note  that  the  isoparametric  formulation  provides  its  elements  with 
the  quality  of  compatibility  and  the  requirements  of  monotonic  convergence. 

The  compatibility  condition  requires  that  the  displacement  along  an  edge  between  the 
adjacent  elements  be  uniquely  defined  in  terms  of  the  displacements  of  the  nodal  points 
along  that  boundary  line  at  any  stage  of  loading  (for  example,  initial  and  final).  Eq.  (3.5) 
and  (3.7)  define  exactly  the  displacements  and  the  coordinates  along  the  edge  of  an  element 
as  functions  of  the  nodal  displacements  and  the  nodal  coordinates,  respectively;  thus 
compatibility  is  ensured. 

A  general  displacement  field  must  allow  an  element  to  undergo  a  rigid-body  displace¬ 
ment  without  introducing  strain.  The  constant  terms  in  the  expression  for  displacement  (for 
example,  Eq.  (3.1))  satisfy  the  rigid- body  translation  while  the  linear  terms  satisfy  the 
rigid-body  rotation. 

A  state  of  constant  strain  will  be  represented  in  an  element  as  its  size  decreases. 

This  condition  is  again  satisfied  by  having  the  complete  linear  expansion  in  the  displace¬ 
ments  (Timoshenko,  1934).  The  premise  is  proved. 

The  property  of  interelement  continuity  and  the  assurance  of  convergence  of  the 
solution  are  important  considerations  for  element  selection.  In  practice,  good  element  per¬ 
formance  invariably  requires  the  exercise  of  sound  judgment  in  locating  a  suitable  set  of 
nodal  coordinates  to  describe  the  geometry  of  the  individual  elements  when  a  structure  is 
idealized.  It  is  important  to  keep  the  element  aspect  ratio  (i.e.,  the  ratio  of  adjacent  edges 
when  an  element  is  proportioned)  from  becoming  excessive.  Also  it  is  essential  to  locate  the 
side  nodes  of  a  curved  edge  close  to  the  center  of  the  edge  and  to  form  corner  angles  well 
under  180  deg.  The  rationale  will  become  clear  through  an  evaluation  of  the  computational 
process.  For  instance,  the  numerical  integration  of  an  isoparametric  element  will  involve 
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calculation  of  the  Jacobian  of  the  coordinate  transformation,  or  the  functional  determinant 
(9  )  y }  3) 

|</|  =  |  - - — - —  |  .  To  ensure  that  the  mapping  is  unique  and  that  there  is  a  one  to  one 

<5  (£,  T],  () 

correspondence  of  coordinates  (%,  y,  2)  and  (<f,  77,  £),  it  is  necessary  that  the  Jacobian  de¬ 
terminant  does  not  change  sign  in  the  element  space  (which  implies  \J\  4  0). 

Now  note  that 


\J\  = 


dx 

dy 

d  z 

di 

d£ 

dx 

dy 

d  2 

d  77 

dr/ 

drj 

dx 

dy 

d  2 

d£ 

(3.25) 


where  Pg,  P  ,  P ^  are  vectors  tangent  to  local  coordinate  (f,  77,  C)  lines  (Fig.  3.4).  From  the 
standard  expression 


d  (Vol)  =  dx  dy  d2 

=  \J\  d{dvdt 


(3.26) 


hence  |«/|  =  lim 


A  Vol  xy  2 


Therefore,  the  Jacobian  provides  a  quantitative  measure  by 


A  Vol  £77  C 

which  the  admissibility  of  an  element  geometry  can  be  evaluated.  This  is  a  useful  guide  in 
selecting  the  appropriate  structural  idealization  for  element  application.  (Since  |</|  =  P ^  ■ 

P ^  x  Pg,  the  Jacobian  is  viewed  as  a  scalar  triple  product,  and  numerically  it  equals  the 
volume  of  a  parallelopiped  with  opposing  sides  parallel  to  the  vector  triplet  Pg,  P  ,  P^.) 


3.3  Specialization 

3.3.1  Load  Matrix  for  a  Prescribed  Pressure 

It  is  often  a  tedious  task  to  calculate  the  load  on  the  surface  of  a  complex  structure 
due  to  a  distributed  pressure.  The  difficulty  arises  from  the  lack  of  a  simple  practical  means 
to  describe  a  design  surface.  However,  with  the  evolution  associated  with  the  isoparametric 
curved  element  formulation,  such  a  surface  can  be  numerically  defined  and  load  calculation 
evaluated. 

In  a  finite  element  system,  the  layout  of  the  element  mesh  pattern  that  represents  a 
structure  depends  on  the  manner  in  which  the  loads  are  to  be  carried.  As  a  rule,  loads  are 
prescribed  only  at  the  nodal  points  and  in  the  directions  corresponding  to  displacement  com¬ 
ponents  defined  in  the  global  coordinates  X,  Y,  and  Z.  Sometimes  “statically  equivalent” 
loads  are  used  as  an  expedient  computation.  For  correct  solutions  especially  when  the 
details  of  local  stress  distribution  are  desired,  the  load  calculations  should  follow  the  pro¬ 
cedure  outlined  in  Section  2.2.1.  An  example  is  given  here. 
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ELEMENTARY  AREA 


Figure  3.4 


n  =  SURFACE  NORMAL 

Curved  element  representation  of  a  complex  surface 


Consider  the  general  case  of  a  fluid  pressure  distribution 


P  =  -  q(x,  y,  2)  ■  ri  (3.27) 

where  q{ x,  y,  s)  is  a  scalar  factor  that  describes  the  variation  in  spatial  pressure  and  n  is  the 
external  surface  normal.  In  case  of  hydrostatic  pressure,  q(%,  y ,  2)  reduces  to  q(z).  From 
Eq.  (2.12),  the  equivalent  loading  vector  at  the  ith  node  is 


1F(/)U 


(3.28) 


where  N  (I)  is  the  shape  function  for  node  i  and  ds  is  the  differential  area  of  a  given  surface 
region  ( D )  for  which  the  pressure  has  been  specified. 

Since  the  position  of  any  point  in  an  element  body  is  defined  by  the  coordinate  equation 
(Eq.  (3.5)),  points  on  any  surface  area  can  be  readily  obtained.  For  instance,  by  setting 
£  =  +  1,  we  obtain  the  surface  equations  in  a  parametric  form  for  the  top  and  the  bottom  faces, 
respectively. 

X  =  '2.NiXi  =  X(£,  r,) 

i 

Y  =  SNiYi  -y(£ij)  (3.29) 

i 

Z=lNiZi  =  Z({,  r,) 

i 


The  summation  extends  to  all  nodes  on  a  given  surface.  Nodal  numbers  for  a  given  (top) 
face  and  corresponding  shape  functions  are  shown  in  the  following  nodal  number  labeling  scheme. 


Quadratic  Element 


Cubic  Element 
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For  the  quadratic  element: 

Corner  nodes  1,  2,  3,  4  £.  =  -  1,  i?.  =  -  1. 

=  J  (1  +  ^  €)  (1  +  r,i  t j)  (<f.  £+  rJir)-  l) 
Midside  nodes  5,  6  £■  =  -  1,  n.  =  0. 

’  +  ’  '  i 

Ni-  + 

Midside  nodes  7,  8  A  =  0,  r?.  =  -  1. 
w-=  -^(l-f2)  (l  +  ij.-i?) 

For  the  cubic  element: 

Corner  nodes  1,  2,  3,  4  <f£.  =  -  1,  ??.  =  -  1. 

Ni  =  ^  (l  +  e^)a  +  r?.^)[9(f2  +  ^)-10] 

Side  nodes  5,  6,  7,  8  n.  =  -  —  ,  =  -  1. 

’  ’  ’  'i  +  3  1  + 

*,•  =  ^  (l  +  ^^)(l-,2)(l  +  9,.,) 

Side  nodes  9,  10,  11,  12  £,  =  x  — ,  n.  =  -  1 

9 

V  ^  (1-f2)  (1  +  9^0  (!  +  >?, -9) 


(3.30) 


(3.31) 


With  the  surface  defined  in  curvilinear  coordinates,  i.e.,  S(g,  77)  in  Eq.  (3.29),  its 
normal  vector  can  be  expressed  by 

»  =  <  +  «y  ?  +  (3.32) 

I  A  A  A  | 

e  j  k 


where  Pg  and  P  are  vectors  tangent  to  the  surface  and  in  the  directions  of  the  tj-  and  rf 
coordinate  lines,  respectively  (Fig.  3.4).  An  elementary  area  of  the  surface  is  given  by 


AS  =  Pg  x  P^  Ar, 


The  surface  area  of  a  region  D  is 

s  -  1  |Pf  *  *4  Af  A’ 

A  77  -»  e 

The  area  of  one  complete  face  of  an  element  can  be  found  conveniently  by  a  numerical  integra¬ 
tion  process;  we  have 


S 


X1  f1  1«  |  d£dr, 


NPT  NPT 

2  2 

i- 1  /=  1 


|S| 


ij=  1,2, 


(3.33) 

NPT 


For  a  relatively  smooth  surface,  such  as  the  surface  of  a  cylinder  or  a  skewed  propeller  blade, 
a  three-point  integration  rule  gives  adequate  accuracy  (i.e.,  for  NPT  =  3,  the  error  range  is 
0.03  to  0.12  percent).  Hence  the  numerical  integration  provides  an  effective  way  to  compute 
a  complex  surface  area  as  well  as  its  projections  and  other  surface  characteristics. 

Now  the  nodal  load  components  for  any  node  (/)  due  to  a  pressure  loading  on  a  face  of 
the  element  can  be  expressed  as 


f(/)=/J/J  N(I)  \p\  |»|  d£  dr,  (3.34) 

Expressed  in  quadrature  format, 

=  ij-N(l).nx.Hitf).H.(r,) 

1  j  / 

y0=^?...«(/).»r  ■  Hi(t)  ■  Hjirj) 

The  distributed  load  q..  must  be  evaluated  in  a  pointwise  manner  at  each  integration  point 
i,  j  =  1,  2,  ...  ,  NPT.  As  before,  these  are  the  numbers  of  Gaussian  integration  points  along 
the  f-  and  rj -coordinate  lines,  respectively.  The  final  load  vector  \p]  is  obtained  by  summing 
up  individual  contributions  from  all  elements  attached  to  these  nodes  or  the  network  of  nodes 
(see  Section  2.2.2). 


j  (3.35) 
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3.3.2  Stresses  on  an  Arbitrary  Surface 


For  obvious  reasons,  it  is  conventional  to  give  experimental  stress  data  along  certain 
well-defined  reference  frames  tangent  to  the  surface  area  of  interest.  Expression  of  the  stress¬ 
es  in  a  global  coordinate  system  generally  does  not  give  a  clear  picture  of  the  surface  situation 
for  a  structure  of  arbitrary  shape.  The  stress  calculation  should  therefore  be  put  into  a  format 
that  allows  numerical  results  to  be  readily  assessed,  interpreted,  and/or  compared  to  experi¬ 
mental  values  or  other  known  results. 

In  the  process  of  computing  stiffness  coefficient  matrix  [A7],  values  of  strain  and  stress 
are  computed  at  each  Gaussian  quadrature  point.  These  values  are  expressed  in  terms  of  nodal 
displacement  parameters  (it v ^  «>.),  with  reference  to  the  global  coordinate  ( x ,  y,  z)  system. 

As  evaluated,  these  strain  values  are  stored  on  tape.  Once  the  equations  of  the  structure  sys¬ 
tem  are  solved  and  nodal  displacements  become  known,  stresses  and  strains  can  be  obtained  by 
direct  substitution  into  Eq.  (2.4)  and  (2.6).  For  design  purposes,  these  results  are  sufficient 
to  describe  the  response  of  an  elastic  body  in  many  practical  applications. 

For  structures  of  complex  shapes,  stresses  at  points  other  than  those  quardrature  points 
may  be  desirable.  This  will  require  additional  calculations  at  the  designated  locations  where 
stress  evaluation  is  sought.  Matrix  [S],  shown  in  Eq.  (3.8),  should  be  used  for  this  purpose. 
The  corresponding  process  outlined  in  Section  3.2.2  must  be  repeated.  Once  again,  there  is  a 
need  to  define  the  surface  orientations  from  which  a  set  of  local  coordinates  can  be  chosen  so 
that  the  computed  stresses  can  be  of  value  for  immediate,  meaningful  interpretations. 

We  begin  with  the  definition  of  unit  normal 


A 


n  = 


n 


(3.36) 


_  —  —  A 

where  n  =  PtxP  is  given  in  Eq.  (3.32).  Now  let  P  be  the  unit  tangent  vector  which  is 

^  '  A  —  —  > 

along  an  ^-coordinate  line  and  P  =  P  /\P  \.  Then 

T  =  P^  x  n  (3.37) 


T  will  be  the  third  unit  vector,  completing  a  right-handed  orthogonal  triad  scribed  in  the  body 

A  A  A 

or  attached  to  a  surface.  The  triad  (T,  P  ,  n)  can  be  considered  as  a  local  reference  frame 
(x ',  y ',  z').  See  Fig.  3.5  and  3.6.  Hence,  the  matrix 


[T,  Py  n]  =  [e15  e2,  e3]  (3.38) 

=  m 

where  e  v  e2,  e 3  are  unit  vectors  in  directions  of  local  rectangular  coordinates  (x y ',  z'). 
[0]  is  also  known  as  the  direction  cosines  matrix.  [T,  P  ,  n\  can  be  expressed  in  terms  of 
traditional  directional  cosine  symbols  (f,  m,  n),  namely 
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Figure  3.5  Rotation  of  reference  frame 


Figure  3.6  Element  nodal  incidence  and  local  reference  frame 


(3.39) 


A  A  A  A 

T  =(.1i  +  m1j  +  n1k 

A  A  A  A 

=  £2  i  +  m2  j  +  n2  k 

a  „  a  A  ? 

ft  =  i3  i  +  m3  3  +  n2  k  ’ 


Let  (w  ',  v',  w')  be  the  displacement  components  along  the  local  cartesian  coordinates 
(x',  y  ',  s').  They  can  be  obtained  from  displacements  (u,  v,  w),  which  are  referenced  to  the 
global  cartesian  coordinates,  by  a  single  transformation: 


u ' 

'  u  ’ 

V  ' 

.  =[/?]. 

V 

w 

w 

k  1 

where  [/?]  is  the  rotation  matrix 


[/?]  = 


mj  ftj 

f-2  m2  n2 

f3  m3  ft3 


=  [0]T 


(3.40) 


Following  the  standard  rule  of  differentiation  (see  Eq.  (3.10))  displacement  derivatives  can  be 
written  with  respect  to  local  coordinates  (x ',  y\  s'): 


d  u 

d  v 

dw 

dx 

dy 

d  z 

d  u 

d  v 

dw 

dx' 

dx' 

dx' 

dx' 

dx' 

d  x' 

d  x 

d  x 

dx 

du 

d  v 

d  w 

d  x 

dy 

d  z 

du 

dv 

d  w 

V 

dy' 

~dy' 

dy' 

dy' 

Ty' 

dy 

dy 

dy 

du 

d  V 

d  w 

d  x 

dy_ 

d  z 

du 

d  v 

d  w 

da' 

\ 

1  ^ 

1 

da' 

da' 

da' 

da' 

d  z 

d  z 

d  s 

=  [R]  < 


[u,  V,  w] 


d 

d  z 


(3.41) 
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By  substituting  local  displacements  («',  v ',  w ')  in  place  of  ( u ,  v,  w )  in  Eq.  (3.41),  we  can 
arrive  at  a  convenient  form  of  the  expression  to  calculate  local  strain  components  from  strains 
given  in  global  coordinates: 

dv' 
dx' 

dv ' 

dy' 

dv' 
dz' 

dv 
dx 

d  v 
dy 

dv 
dz 

Stress  components  computed  in  the  global  coordinate  system  can  be  transformed  to 
locally  oriented  stresses  [ct']  by  a  similar  expression,  that  is, 


=  [ff][a][tf]r  (3.43) 

Once  the  element  mesh  and  labeled  element  nodal  numbers  are  laid  out  over  the  ideal¬ 
ization  of  a  structure,  the  element  incidences  must  be  suitably  ordered  (Fig.  3.6)  such  that 
one  curvilinear  coordinate,  77  for  instance,  will  be  placed  on  a  designated  coordinate  surface. 
When  the  surface  normal  vector  is  computed,  another  orthogonal  surface  tangent  will  complete 
a  right-handed  triad.  This  combination  will  furnish  a  set  of  local  rectangular  axes  and  form 
the  basis  of  a  rotation  matrix.  These  data  enable  the  stresses  to  be  calculated  over  an 
arbitrarily  shaped  body  along  any  prescribed  orientation . 
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3.3.3  Applications  to  Plates  and  Shells 


The  solid  isoparametric  finite  elements  outlined  here  are  based  on  a  general  three- 
dimensional  solution.  Such  elements  have  demonstrated  broad  applicability  for  problems  in 
structural  mechanics.  Where  the  element  thickness  is  decreased  to  the  proportion  of  a  med¬ 
ium  thick  plate  (or  a  thin  shell),  a  specialized  formulation  is  permissible  to  achieve  greater 
economy  and  effectiveness  (Ahmad  et  al.,  1965,  1970  and  Pawsey,  1970).  Some  well-known 
approximation  will  be  utilized  in  the  computations  of  stiffness  coefficients  [A'  ],  e.g.,  as  ad¬ 
vocated  in  classical  plate  theory  (Timoshenko  and  Woinowsky-Krieger,  1959  and  Flugge, 
1960)  such  that  lines  perpendicular  to  the  middle  surface  remain  straight  under  loading  and 
strains  along  these  lines  can  be  ignored  in  the  energy  summation. 

We  begin  with  the  element  geometry.  Pairs  of  points  £top  and  *bottom,  each  with  given 
cartesian  coordinates,  prescribe  the  shape  of  an  element  (see  Fig.  3.7)  i=  1,  2,  ...  ,  NNPE 
where  NNPE  is  the  number  of  nodes  per  element.  Let  f,  rj  be  the  two  curvilinear  coordinates 
in  the  middle  surface  of  the  shell  and  £  a  linear  coordinate  in  the  thickness  direction.  As 
before,  77,  and  £  vary  between  - 1  and  +  1  on  the  respective  faces  of  the  element.  The 
cartesian  coordinates  (x,  y ,  z)  of  any  point  of  the  shell  can  be  defined  by  interpolation  on 
the  coordinates  at  nodal  points  i\  hence 


x 

y 

z 


X  N. 

i'=i  1 


2  N. 
1=1  1 


(1-0 

2 


Vi 

Z  ; 


bottom 


(3.44) 


Here  /V.  =  AL(f,  77)  are  surface  shape  functions  of  the  type  given  by  Eq.  (3.30).  Once  again, 
parabolic,  cubic,  etc.  (or  of  any  specific  order)  shape  functions  can  be  adopted  for  the  middle 
surface  of  a  shell  element. 

By  introducing  a  nodal  vector  V  that  connects  the  pairs  of  nodes  i.  and  i.  ..  ,  we 

1  top  bottom’ 

can  rewrite  the  relationship  between  the  cartesian  and  curvilinear  coordinates  in  terms  of  the 
midsurface  coordinates  and  the  vector  V, .. 


(3.45) 
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where 


^  =  < 


X  . 

I 


Vi 


si 


x . 

i 


Vi 


3i 


(3.46) 


top  \  /  bottom 

<=  1,2,  ,  NNPE  (=n) 


Now  the  displacement  pattern  has  to  be  assumed  for  the  element.  Since  the  strains  in 
the  direction  normal  to  the  midsurface  are  assumed  to  be  negligible,  the  displacement  vector 
r  throughout  the  element  will  be  taken  to  be  completely  defined  by  the  midsurface  nodal  dis¬ 
placement  r  ■  and  two  rotations  (a,  /3)  of  the  nodal  normal  V, .  about  axes  orthogonal  to  it.  If 

A  A 

two  such  orthogonal  directions  are  given  by  unit  vectors  V2[.  and  Vx with  corresponding 
scalor  rotations  a  -  and  B.,  we  can  write 


where 


and 


N. 


I 


(a.  +  /3.)  x  V 


3  i 


^  be¬ 
lli  +  V  ]+  w  k 


w 


T  .  -  U.  I 
m  i  i 


) 


A  A 

+  w  ■  k 


(3.47) 


(3.48) 


! 


A 


I  (3.50) 


*=  1,2,  ...  ,  NNPE  (=  n) 


Here  u,  v,  and  w  are  displacements  in  the  directions  of  global  x,  y,  and  z  axes,  and  u.,  v and 
w i  are  displacements  at  the  midsurface  node  i.  At  each  node,  we  now  have  the  five  basic 
degrees  of  freedom: 
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(3.51) 


From  the  basic  shell  assumptions,  the  strain  components  are  essential  in  directions 
of  tangential  orthogonal  axes  related  to  the  surface  C  =  constant.  The  strains  to  be  used  in 
the  element  approach  must  be  converted  to  this  same  reference  system.  At  a  point  on  this 
surface,  take  z '  in  the  direction  of  surface  normal  n,  Eq.  (3.32).  We  can  establish  a  set  of 
local  orthogonal  axes  x ',  y',  and  z One  simple  scheme  is  given  by  Eq.  (3.49).  The  strain 
components  of  interest  are 


■y  / 

7  x  y 


y  ✓  / 
'  x  z 


du' 

dx' 

dv' 

dy' 

du' 

du' 

~dy' 

+ 

~dx' 

du' 

dw' 

~dl' 

+ 

lx' 

dv' 

+ 

dw' 

Jz' 

ly' 

[B']  W“\ 


(3.52) 


The  stresses  corresponding  to  these  strains  are  defined  with  the  aid  of  the  elasticity 
matrix  [D  '].  We  have 


=  W']  ir 


(3.53) 
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where 
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(3.54) 


The  5x5  matrix  [D  ']  is  defined  for  an  isotropic  material.  E  and  v  are  Young’s  modu¬ 
lus  and  Poisson’s  ratio,  respectively.  The  factor  k  (=  1.2)  included  in  the  last  two  shear 
terms  is  intended  to  improve  the  shear  displacement  approximation  (Ahmad,  1970).  It  is  seen 
that  because  of  the  displacement  assumption,  the  shear  strain  is  approximately  constant 
through  the  thickness,  whereas  in  reality  the  shear  distribution  is  approximately  parabolic. 
The  value  k=  1.2  is  the  ratio  of  relevant  strain  energies. 

The  next  step  is  the  calculation  of  the  element  stiffness  matrix 


[K]  =  f1  f1  f1  [B']t[D'\  [B']  \J\  dr,  d(  (3.55) 

—l—l  — l 

By  definition,  =  [S']  Matrix  [S']  relates  the  local  displacement  derivatives  to  the 

nodal  parameters.  The  calculation  of  [B '}  involves  three  steps: 

a.  Compute  global  displacement  derivatives  for  a  set  of  curvilinear  coordinates  in  the 
manner  shown  by  Eq.  (3.10). 

b.  Transform  these  derivatives  (i.e.,  the  strain  components  expressed  in  global  coordi¬ 
nates)  to  local  displacement  derivatives  by  Ea.  (3.42).  Here,  the  direction  cosine  matrix 

A 

[0]  can  be  constructed  by  a  process  given  by  Eq.  (3.49)  with  unit  vector  F3  parallel  to  z  -axis 
which  is  in  the  direction  of  surface  normal  n. 

c.  Assemble  the  local  displacement  derivatives  to  form  the  strain  vector  [r'l  in  terms  of 
nodal  parameter  vector  ft/°|  given  by  Eq.  (3.52). 

Now  the  whole  integral  of  Eq.  (3.55)  can  be  expressed  as  an  explicit  function  of  the 
curvilinear  coordinates.  After  carrying  out  some  operations  at  the  submatrix  level,  simplifi¬ 
cations  and  saving  in  numerical  processing  can  be  achieved.  A  numerical  integration  will 
allow  the  properties  of  the  element  to  be  evaluated. 
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3.4  Implementation 

3.4.1  Introduction  to  Solution  Methods 

In  the  displacement  method  of  finite  element  analysis,  the  problem  eventually  reduces 
to  the  solution  of  a  set  of  linear  simultaneous  equations  that  express  the  load-displacement 
or  equilibrium  relation  for  the  structure.  The  displacement  boundary  conditions  can  be  read¬ 
ily  imposed  by  deleting  the  appropriate  nodal  parameters  (corresponding  to  nodal  degrees  of 
freedom)  from  Eq.  (2.14).  After  reduction— which,  in  effect,  removes  the  rigid-body  mode— we 
have 


[Kr]  V  =  P  (3.56) 

where  [K r]  is  the  nonsingular  structure-stiffness  matrix, 

U  is  a  vector  of  unknown  nodal  parameters  (displacements),  and 

JP  is  a  vector  of  applied  load. 

Once  Eq.  (3.56)  is  solved,  the  global  displacement  parameters  U  become  known.  The 
desired  stress  and  strain  at  any  point  within  any  element  can  be  found  immediately  by  sub¬ 
stituting  the  nodal  displacements  in  Eq.  (2.4)  and  (2.5)  in  turn. 

In  practical  application,  a  sizable  number  of  finite  elements  is  required  for  the  repre¬ 
sentation  of  structural  design  problems.  Consequently,  an  extensive  network  of  nodes  evolves 
and  the  size  of  the  stiffness  matrix  which  corresponds  to  the  number  of  unknown  nodal  vari¬ 
ables  is  often  overwhelming  (not  infrequently  several  thousand  degrees  of  freedom  arise).  A 
large  portion  of  the  total  computer  time  required  to  solve  a  given  problem  is  generally  con¬ 
sumed  in  solving  the  set  of  linear  equations  (Eq.  (3.56)).  Here,  the  method  of  solution  can 
have  a  significant  bearing  on  the  computational  efficiency  which  is  measured  in  terms  of 
demand  on  core  size.  At  times,  the  core  requirement  may  dictate  the  applicability  of  the 
finite  element  method. 

It  is  of  prime  interest,  then,  to  select  a  solution  algorithm  which  takes  into  account 
the  symmetric,  positive  definite  and  banded  nature  of  stiffness  matrix  [K  ].  Further  the 
Gauss  elimination  is  known  to  be  numerically  stable,  irrespective  of  the  order  in  which  the 
equations,  Eq.  (3.56),  are  eliminated  (pivot  search  is  not  necessary)  and  therefore  the  full 
advantage  of  symmetry  can  be  realized.  The  elimination  of  a  row  S  (which  represents  an 
equilibrium  equation  in  nodal  variable  U$)  leads  to  a  modification  of  the  coefficients  in  the 
remaining  rows  according  to  the  formulas 


(3.57) 
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Stiffness  coefficients  k..,  k.  (=  k  .),  .  .  etc.  represent  the  sum  of  individual  element  contri- 
butions.  It  does  not  matter  in  which  order  the  summation  is  made  and,  further,  they  need  not 
be  fully  summed  except  those  in  row  S  currently  being  eliminated.  This  process  affects  a 
triangular  array  immediately  below  the  row  that  is  being  eliminated.  It  can  be  carried  out  by 
retaining  in  the  core  only  the  triangle  of  coefficients  which  moves  diagonally  downward  as 
the  elimination  proceeds;  see  Fig.  3.8. 

As  depicted  by  a  coarse  or  a  fine  mesh,  the  finite  element  idealization  of  a  structure 
frequently  takes  the  form  of  a  simply  connected  region.  In  these  cases,  the  stiffness  matrix 
can  always  be  arranged  in  a  nicely  banded  form  by  labeling  the  unknown  parameters  U ■  in  a 
suitable  order.  In  the  case  of  a  closed  ring  that  represents  a  multiply  connected  region,  the 
band  can  be  large  despite  high  sparsity.  It  is  also  true  that  the  bandwidth  of  the  individual 
equations  tends  to  be  large  in  three-dimensional  problems  where  large  solid  elements  with 
intermediate  nodes  along  their  edges  are  used.  In  addition,  one  has  to  operate  on  many  zero 
terms  within  the  band,  and  this  adds  to  the  cost  of  solving  equations  by  a  band  algorithm. 

For  these  and  many  other  cases  as  well,  the  “frontal”  technique  (see  Section  3.42)  requires 
less  storage  and  arithmetic  than  the  best  Gaussian  banded  algorithm  (Melosh  and  Bamford, 

1969  and  Irons,  1970).  The  frontal  solution  algorithm  is  employed  for  this  study,  and  is 
described  in  the  next  section. 

3.4.2  Frontal  Technique 

Most  of  the  inefficiency  in  data  processing  keyed  to  matrix  bandedness  can  be  avoided 
by  applying  the  frontal  technique.  It  is  based  on  the  principle  implied  by  the  very  nature  of 
Gauss  elimination,  Eq.  (3.57).  Frontal  processing  takes  advantage  of  the  fact  that  nonzero 
elements  in  a  column  of  the  decomposition  (inside  the  triangular  portion  of  the  array  affected 
by  row  operations)  cannot  occur  in  any  row  prior  to  the  occurrence  of  a  nonzero  element  (i.e., 
the  stiffness  coefficient  kis  =  kg.  =/  0)  in  the  column  of  the  stiffness  matrix.  In  this  approach, 
the  first  appearance  during  decomposition  of  an  element  in  a  column  of  [K r ]  causes  the  addition 
of  that  column  to  the  wavefront  which  is  the  “front  of  active  variables.”  Data  that  must  be 
readily  available  on  core  include  only  the  coefficients  of  the  equations  on  the  wavefront.  A 
variable  becomes  “active”  on  its  first  appearance  and  is  eliminated  immediately  after  its 
last.  Hence,  if  a  variable  xg  is  ready  for  elimination,  there  must  be  no  subsequent  elements 
that  contain  x  .  (It  was  implied  previously  that  kg.,  A  .,  .  .  .  .  ,  Pg  in  Eq.  (3.57)  must  be 
fully  summed.)  Thus  the  size  of  the  “front”  is  smaller  than  or,  at  worst,  equal  to  the  band¬ 
width  of  the  equations.  In  effect,  this  eliminates  a  good  number  of  zeroes  from  the  process 
of  computation,  and  it  saves  substantially  on  core  requirement. 

Fig.  3.9  illustrates  an  application  of  the  wavefront  method.  The  associated  graph, 
representing  a  portion  of  a  ship  structure,  explains  the  name  “wavefront,”  i.e.,  equation 
solving  is  visualized  as  progressing  like  a  wave  over  the  structure.  At  any  given  time  in 
the  analysis,  the  wavefront  includes  the  total  number  of  active  degrees  of  freedom. 
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FRONT _ _ 

PROGRESSING  NUMBER  OF  ELEMENTS  =  17 

NUMBER  OF  JOINTS  =  85 


LEGEND  BOUNDARY  NODES 

(M)  =  ELEMENT  NUMBER  ®  Ui  =  0 

N  =  JOINT  NUMBER  @  Uj=Vj  =  0 

Figure  3.9  Frontal  processing  of  a  finite  element  idealization  of  the  cross  frame  of  a  ship 


50 


These  are  associated  with  partially  processed  degrees  of  freedom  (or  nodal  variables).  For 
example,  Fig.  3.9  shows  a  situation  just  before  element  (8)  is  introduced.*  If  there  is  only 
one  equation  per  joint,  joints  7,  17,  27,  and  60,  70,  77  (representing  active  nodal  variables) 
currently  constitute  the  wavefront.  Joints  1,  2,  3,  4,  5,  6,  14,  15,  16,  21,  22,  23,  24,  and  25, 
etc.  (which  lie  to  the  left  (back)  of  the  dotted  line)  have  been  fully  processed  and  are  no 
longer  considered  active.  At  the  next  step,  element  (8)  introduces  the  new  active  variables 
61,  62,  71,  78,  and  79.  As  elements  are  being  processed  one  after  another,  the  front  moves 
forward  across  the  structure. 

Further  we  observe  that: 

a.  The  order  of  element  sequence  is  critical  in  frontal  solution,  just  as  the  node  numbering 
is  critical  in  a  band  algorithm,  although  the  best  order  is  not  always  unique. 

b.  Ordering  of  the  nodes  is  irrelevant  to  the  frontal  technique.  The  node  numbers  are 
merely  unique  labels  that  relate  the  degrees  of  freedom  for  element  and  structure.  They  have 
no  effect  on  the  order  of  elimination.  Here,  when  a  given  mesh  is  changed  by  adding  or  de¬ 
leting  nodes,  the  frontal  data  are  little  changed,  but  a  band  algorithm  may  require  extensive 
node  renumbering  in  order  to  preserve  a  small  bandwidth.  This  can  be  extremely  valuable  in 
design  application  when  a  mesh  is  to  be  modified  locally  for  a  rerun. 

c.  The  size  of  front  and  hence  the  storage  requirements  for  a  given  problem  can  be 
assessed  from  the  mesh  pattern  such  as  given  in  Fig.  3.9.  For  example,  before  element 

(8)  is  introduced,  there  are  six  active  variables:  7,  17,  27,  60,  70,  77;  immediately  after 
the  stiffness  coefficient  generated  by  element  (8)  has  been  added  in,  the  list  of 
active  variables  increases  to  11. 

3.5  Evaluation  of  Numerical  Results 

The  preceding  sections  have  described  the  formulation  of  solid  finite  elements  for  the 
three-dimensional  analysis  of  a  general  elastic  body  and  a  method  for  its  implementation. 

The  performance  of  these  elements  in  the  representation  of  an  elastic  structure  will  now  be 
evaluated.  Some  preliminary  computations  were  conducted  with  these  elements  using 
examples  given  in  Chapter  9  of  Zienkiewicz  (1970)  and  favorable  numerical  results  were 
obtained.  Additional  problems  are  included  here  and  solutions  for  such  basic  structural 
components  as  beams,  plates,  and  shells  are  presented  to  verify  the  behavior  of  individual 
elements  and  to  examine  the  adequacy  of  the  finite  element  models.  A  special  example  will 
be  given  in  Chapter  4  to  demonstrate  the  effectiveness  of  the  procedure  in  solving  the 
difficult  problem  of  marine  propeller  blades. 


*  Element 


has  not  yet  been  introduced. 
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3.5.1  Prismatic  Beams 


3.5. 1.1  A  Cantilever  Beam.  Timoshenko  (1934)  treated  a  prismatic  bar  subjected  to  pure 
bending  as  a  problem  of  elasticity  in  three  dimensions.  Consider  the  bar  bent  in  one  of  its 
principal  planes  (x-z)  by  a  couple  M  at  one  end  and  held  fixed  at  the  other.  Take  the  origin 
of  the  coordinates  at  the  centroid  of  the  cross  section,  i.e.,  point  0  (Fig.  3.10).  The  displace¬ 
ment  components  as  predicted  by  the  theoretical  solution  are: 
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2  R 


[z2  -  v  (x2  -  y2)] 
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1  M 

and  • —  =  —  is  the  curvature  of  the  deflected  bar.  After  deformation,  points  on  a  plane 
R  El 

cross  section  will  remain  on  that  plane,  but  the  cross-sectional  shape  will  change.  The 
sides  of  the  cross  section  become  inclined  and  the  top  and  bottom  edges  are  bent  into 
parabolic  curves  as  shown  by  the  dotted  lines  in  Section  A- A  of  Fig.  3.10. 

An  aluminum  cantilever  beam  having  a  cross  section  of  0.5-  by  0.6875-in.  thick  and 
a  length  of  11.0  in.  was  bent  by  an  end  couple  of  4.10  in. -lb.  Surface  displacements  were 
measured  by  using  holographic  interferometry  (Dhir  and  Sikora,  1972).  Three-dimensional 
displacement  components  («,  v,  w)  for  the  top  face  along  the  front  edge  were  obtained;  see 
Fig.  3.10.  This  experiment  was  conducted  as  a  preliminary  process  to  evaluate  the  quality 
of  measurements  obtainable  by  holography.  In  this  case,  the  measured  data  indicate  that  the 
beam  appeared  to  be  distorted  locally  near  the  support;  the  investigators  had  difficulty  in 
simulating  the  fixed-end  condition  for  the  slender  beam. 

The  cantilever  beam  is  represented  by  four  20-node  (quadric)  hexahedron  elements 
with  a  total  length  of  11  in.  The  finite  element  model  is  completely  fixed  at  one  end  and 
bent  at  the  other  by  an  end  couple  produced  by  two  equal  and  opposite  forces  (=+  5.9636  lb) 
applied  to  the  top  and  bottom  midside  node.  Computed  displacements  u  and  v  are  almost 
exact;  they  are  within  0.1  percent  of  the  correct  solution.  Computed  stresses  are  good  and 
are  generally  within  0.5  percent  of  the  exact  solution  except  at  points  near  the  concentrated 
applied  forces. 


3.5. 1,2  A  Simply  Supported  Beam.  The  purpose  of  the  current  example  is  to  assess  the  effect 
of  mesh  pattern  (as  it  relates  to  the  element  aspect  ratio)  on  the  behavior  of  the  finite  ele¬ 
ment  model.  A  simply  supported  rectangular  beam  40  in.  long  and  with  a  0.6  by  6.0  iri.  cross 
section  is  represented  by  four  quadric  elements  (Table  3.1). 
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EMENT  (INCH  X  103) 


TABLE  3.1  RESULTS  OF  ANALYSIS  OF  A  SIMPLY  SUPPORTED  BEAM 


DISCRETE  MODEL  B1 


DISCRETE  MODEL  B2 


SPAN  £  =  40  IN. 

THICKNESS  t  =  6  IN. 

WIDTH  b  =  0.6  IN. 

CENTER  POINT  LOAD  P  =  10,000  LB 


POISSON  RATIO  V  =  0.30 

YOUNG'S  MODULUS  E  =  30,000,000  PSI 

BENDING  STRESS  a  PSI 


Vertical  Deflection 

at  Section 

in. 

Stresses  at 

Section  C 

psi 

A 

B 

C 

D 

°top 

^bottom 

Discrete  Model  B1 

0.0156 

0.0292 

0.0390 

0.0427 

-  20,590 

21,030 

Discrete  Model  B2 

0.0155 

0.0290 

0.0388 

0.0425 

-  20,820 

21,040 

Beam  Theory 

0.0151 

0.0283 

0.0376 

0.0412 

-  20,830 

20,830 

Theory  of  Elasticity 

0.0158 

0.0298 

0.0398 

0.0440 

-20,720 

20,890 
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These  elements  are  of  uniform  size  in  discrete  model  B1  while  they  are  of  nonuniform 
size  in  discrete  model  B2.  The  beam  is  loaded  at  midspan  by  a  concentrated  force  P . 

Vertical  deflections  and  normal  stresses  are  obtained  from  the  discrete  element  models.  The 
response  of  discrete  model  B2  with  its  nonuniform  mesh  size  appears  to  be  slightly  stiffer. 
Theoretical  solutions  (listed  in  Table  3.1)  are  obtained  (a)  by  using  elementary  beam  theory 
in  which  the  effect  of  shear  stresses  is  neglected  and  (b)  by  elastic  theory  (Timoshenko  1934) 
in  which  the  beam  is  treated  as  a  two-dimensional  problem.  Numerical  results  obtained  by 
both  finite  element  calculations  are  acceptable  from  an  engineering  standpoint. 

3.5.2  Plate  Bending 

When  the  thickness  of  an  elastic  body  is  small  compared  to  its  other  dimensions,  it 
is  called  a  plate.  A  flat  plate  is  classified  as  thin  or  medium  thick  according  to  its  span 
to  thickness  ratio  (t./t)  whether,  this  ratio  is  greater  or  less  than  40.  When  a  medium  thick 
plate  (5  >  i/t  >  40)  is  bent  with  small  deflection— small  in  the  sense  that  the  plate  deflec¬ 
tion  is  small  compared  with  its  thickness— the  classical  plate  theory  (i . e. ,  the  Kirchoff 
assumption)  holds.  Then,  the  laterally  loaded  plate  can  be  treated  as  a  two-dimensional 
problem. 

The  20-node  hexahedron  elements  are  used  to  model  a  plate  structure  with  varying  mesh 
fineness  (Fig.  3.11).  Here  the  analysis  is  for  a  simply  supported  square  plate  with  a  40-in. 
span  and  a  1-in.  thickness  subjected  to  a  concentrated  center  load.  This  is  a  severe  bending 
test  because  the  situation  is  singular  at  the  loading  point  (the  solution  for  stresses  at  this 
point  is  not  bounded). 

The  center  deflection  under  concentrated  load  obtained  from  finite  element  calcula¬ 
tions  is  plotted  against  mesh  size  in  Fig.  3.12.  Note  the  rapid  convergence  of  the  deflec¬ 
tion.  For  a  mesh  size  /V  =  4  (t/t  =  5),  the  calculated  deflection  is  within  1.5  percent  of  the 
exact  solution.  Numerical  results  for  some  well-known  plate  bending  elements  are  also 
shown  in  Fig.  3.12  for  comparison  purposes.  These  results  are  taken  from  the  Nastran 
Theoretical  Manual  (1969). 

The  dotted  lines  in  Fig.  3.13  indicate  the  distribution  of  bending  moments  M x  and  M 
along  the  centerline  of  the  square  plate  computed  from  the  finite  element  model.  The  theo¬ 
retical  results  (Timoshenko,  1959)  are  for  a  concentrated  load  P  applied  over  a  circular  area 
of  0.05-a  radius. 

For  general  plate  bending  analysis,  the  application  of  the  specialized  element  given 
in  Section  3.3.3  is  preferred  because  of  its  efficiency  and  ensuing  economy.  Although  the 
solid  element  is  not  intended  to  be  used  for  plate  bending  problems,  it  did  behave  well  as  a 
plate  when  molded  to  the  geometry  of  a  plate.  The  element  did  yield  results  which  converged 
rapidly  to  the  exact  solution. 
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a  =  20  IN.  P  =  4,000  LB 

t  =  1  IN.  E  =  30,000  KSI 

v  =  0.30 

Figure  3.11a  A  simply  supported  square  plate  under  center  load 


Figure  3.11b  Element  mesh  used  on  the  analyzed  quadrant  of  the  plate 


Figure  3.11  Plate  bending  sample  problem  (example  3.1) 
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MOMENT  COEFFICIENT 
P  (IN.-LB/IN.) 


3.5.3  Thick-Walled  Cylinder 

The  subject  of  this  analysis  is  a  section  of  cylindrical  wall  12-in.  long  with  an  outside 
diameter  of  20  in.  The  wall  is  2-in.  thick  and  its  outer  surface  is  under  a  uniformly  distributed 
pressure.  The  problem  is  that  of  an  axisymmetric  solid  and  the  situation  is  two-dimensional. 

In  this  case,  a  quadrant  of  the  thick-walled  cylinder  is  considered  and  each  time  either  two  or 
three  20-node  hexahedron  elements  are  used  to  model  the  curved  shell. 

Displacements  are  computed  for  all  nodal  joints  in  global  coordinates  XYZ.  Stresses 
are  computed  at  integration  points  as  well  as  on  the  curved  surfaces,  including  both  interior 
and  exterior.  The  stresses  are  expressed  in  global  coordinates  and  also  in  local  curvilinear 
coordinates  (£,  77,  £)•  A  two-element  representation  of  the  curved  shell  is  illustrated  in  Table 
3.2.  Also  given  there  are  computed  stresses,  the  radial  component  a,  and  the  tangential  com¬ 
ponent  oq  along  a  typical  section.  Numerical  results  obtained  from  the  three-dimensional 
finite  element  solution  are  in  good  agreement  with  the  theoretical  solution  (Timoshenko,  1934). 
The  computed  stresses  are  generally  closer  to  the  exact  solution  at  the  integration  points  than 
at  other  points,  e.g. ,  on  the  surfaces. 

This  example  provides  a  way  to  verify  certain  subprograms  used  to  define  surface  char¬ 
acteristics,  local  axis  orientations,  and  stress  transformations  for  an  arbitrary  shaped  body. 
These  subprograms  are  designed  to  enable  immediate  interpretation  of  finite  element  solutions 
where  a  complex  structural  configuration  is  considered  (see  Section  4.4). 

3.5.4  Stiffened  Plates 

In  the  case  of  a  simple  beam  having  a  longitudinal  plane  of  symmetry  with  external 
forces  also  acting  in  this  plane,  bending  will  take  place  in  the  same  plane.  Elementary  beam 
theory  assumes  that  a  transverse  section  of  the  beam  that  is  originally  plane  remains  plane 
and  normal  to  the  longitudinal  fibers  of  the  beam  after  bending.  Direct  measurement  shows 
that  the  beam  theory  gives  very  accurate  results  for  the  deflection  of  beams  and  the  strains 
of  longitudinal  fibers.  The  results  of  finite  element  calculation  for  a  beam  of  rectangular 
section  (Table  3.1)  also  confirm  the  validity  of  the  elementary  theory  of  bending. 

Now  consider  a  T-beam  which,  in  essence,  is  a  rectangular  beam  with  projected  flanges. 
When  a  T-beam  has  the  proportions  of  an  ordinary  structural  shape  as  found  in  a  standard  hand¬ 
book,  its  behavior  under  simple  bending  can  still  be  predicted  by  ordinary  beam  theory.  How¬ 
ever,  when  the  flanges  of  a  beam  are  extended  further,  the  transverse  stress  distribution  is  no 
longer  uniform  across  the  width  of  the  flange.  It  is  known  that  parts  of  the  flange  at  a  dis¬ 
tance  from  the  web  do  not  take  their  full  share  in  resisting  bending.  This  phenomenon  is 
called  “shear  lag.”  A  beam  of  this  type  is  less  stiff  than  predicted  by  elementary  beam 
theory. 

The  configuration  of  a  beam  with  very  wide  flanges  becomes  that  of  a  stiffened  plate. 
Plates  that  are  stiffened  transversely  and/or  longitudinally  are  important  components  in  the 
hull  structure  of  ships.  The  plate  which  serves  as  load  bearer  also  acts  as  the  flange  of 
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TABLE  3.2  RESULTS  OF  STRESS  CALCULATIONS  FOR  A  THICK-WALLED  CYLINDER 
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9-000  0.54  4.96  0.58  4.97 


the  stiffener.  There  are  numerous  approximate  methods  for  analyzing  these  plates.  Most  fre¬ 
quently,  they  are  based  on  idealization  of  the  system  into  a  grid  (or  network)  of  beams  or  into 
an  equivalent  orthotropic  plate.  These  methods  require  an  estimate  on  the  effective  width  of 
the  plate  and  the  torsional  rigidity  of  the  composite.  Thus  they  demand  considerable  engi¬ 
neering  judgment  and  interpretation  and  are  not  completely  reliable. 

Important  progress  has  been  made  in  recent  years  on  finite  element  analysis  of  eccen¬ 
trically  stiffened  plates  and  shells  (Mehrain,  1967;  Kohnke,  1969;  and  Chu,  1970).  The  stiff¬ 
ened  plates  are  represented  by  a  composite  of  beams  and  plates  (one-  and  two-dimensional 
elements,  respectively).  Actually,  at  or  near  the  interface  of  the  plate  and  its  stiffener,  an 
accurate  solution  of  their  structural  behavior  would  require  a  three-dimensional  analysis. 

The  improved  hexahedron  element  is  utilized  here  to  analyze  the  intricate  stress  distribution 
of  stiffened  plates. 

Fig.  3.14  indicates  a  general  scheme  used  to  analyze  a  laterally  loaded  stiffened  plate. 
The  plate  with  a  stiffener  on  one  side  is  taken  as  a  simply  supported  plate  beam  carrying  a 
center  load  P.  Because  of  symmetry,  a  half  span  of  the  beam  is  idealized  and  twelve  20-node 
hexahedron  elements  are  used  for  its  representation.  At  each  end,  vertical  supports  are  pro¬ 
vided  for  both  flange  and  web  plates.  Five  different  flange  versus  stiffener  aspect  ratios  are 
used  as  parameters  (plate  beams  CX2  through  CX6)  to  assess  the  effect  of  flange  width 
on  the  structural  behavior  of  these  stiffened  plates.  As  the  flange  of  a  beam  is  progressingly 
widened  (e.g.,  plate  beam  CX3),  the  maximum  center  deflections  calculated  by  the  finite 
element  representation  are  found  to  increase  at  a  greater  rate  than  those  obtained  by  the 
elementary  beam  formula.  In  other  words,  because  of  shear  lag,  each  incremental  flange 
material  is  engaged  in  a  lesser  capacity  than  has  been  assumed. 

Fig.  3.15  depicts  the  distribution  of  longitudinal  stresses  (a  )  on  a  transverse  section 
of  plate  beam  CX3.  Stresses  are  generally  higher  in  the  parts  of  flange  near  the  stiffener 
and  especially  in  the  vicinity  of  the  concentrated  load.  Similar  patterns  of  stress  distribution 
with  varying  degree  of  stress  gradient  are  found  in  each  beam  that  has  projected  flanges. 

Stress  distribution  across  the  web  plate  differs  somewhat  from  the  linear  variation  assumed 
in  the  elementary  theory  of  bending.  The  neutral  axis  at  these  cross  sections  of  the  beam 
does  not  pass  through  the  geometric  centroid  (Fig.  3.15).  Such  behavior  was  more  pronounced 
for  a  plate  beam  such  as  CX3  with  its  deeper  web  plate  (S=  6  in.)  than  for  the  others,  for 
example,  CX4. 

Fig.  3.15  also  shows  the  distribution  of  typical  transverse  stresses  (a)  along  the 
flange  plate.  They  indicate  some  local  bending  of  the  flange  in  the  transverse  direction  near 
the  concentrated  force.  Additional  stress  data  of  interest  are  given  in  Table  3.3.  These 
stresses  were  obtained  at  the  Gaussian  integration  points  (used  to  form  the  stiffness  matrix 
of  the  isoparametric  element).  Stresses  computed  at  such  points  have  been  shown  to  be 
generally  of  high  accuracy.  Stresses  for  other  locations  in  the  plate  beam  can  be  found  by 
interpolation. 
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TABLE  3.3  LONGITUDINAL  STRESSES  (cy)  IN  A  STIFFENED  PLATE  (EXAMPLE  3.2) 


'  Table  3.3a  Plate  beam  CX2 


V 

oy  (psi) 

(in.) 

A1 

A2 

A3 

B1 

B2 

B3 

Cl 

C2 

C3 

0.676 

-6252 

-5103 

-3479 

-  6795 

-5417 

-3525 

-7768 

-6215 

-  4127 

3.000 

-5989 

-4945 

-3529 

-  5845 

-4817 

-3450 

-5959 

-5000 

-  3755 

5.324 

-5557 

-4471 

-3074 

-  5267 

-4442 

-3411 

-5063 

-4555 

-  3961 

9.500 

-3683 

-3180 

-2547 

-  3725 

-3201 

-2537 

-3748 

-3199 

-  2498 

16.500 

-1155 

-  979 

-  752 

-  1214 

-1039 

-  814 

-1255 

-1082 

-  863 

D1 

D2 

D3 

El 

E2 

E3 

FI 

F2 

F3 

0.676 

-9161 

-6841 

-4624 

-10247 

-7926 

-5705 

-2533 

5995 

18193 

3.000 

-6093 

-5051 

-4045 

-  6575 

-5528 

-4519 

-2225 

5486 

16420 

5.324 

-4773 

-4496 

-4048 

-  4646 

-4367 

-3916 

-2135 

4687 

14265 

9.500 

-3757 

-3232 

-2466 

-  3639 

-3113 

-2347 

-1126 

3541 

9634 

16.500 

-1200 

-1050 

-  854 

-  1230 

-1079 

-  882 

-  387 

1203 

3182 

2 H  2 


POSITION  OF  INTEGRATION  POINTS  USED  IN  STRESS  TABULATION 


Tabic  3.3b  Plate  beam  CX3 


Y 

Oy  (psi) 

(in.) 

A1 

A2 

A3 

B1 

B2 

B3 

Cl 

C2 

C3 

0.676 

-3267 

-2335 

sra 

-5092 

-3791 

-  2480 

3.000 

-3428 

-2378 

B39 

-3878 

-2995 

-  2119 

5.324 

-3454 

-2267 

-1076 

-3158 

-2311 

-1473 

-3124 

-2641 

-  2177 

9.500 

-2015 

-1629 

-1244 

-2151 

-1717 

-1283 

-2268 

-1803 

-  1338 

16.500 

-  522 

-  369 

-  218 

-  729 

-  586 

-  444 

-  805 

-  670 

-  537 

D1 

D2 

D3 

El 

E2 

E3 

FI 

F2 

F3 

0.676 

-6733 

-4610 

-2802 

-7851 

-5731 

-3925 

-1290 

+5261 

14095 

3.000 

-4173 

-3260 

-2491 

-4656 

-3745 

-2980 

-1004 

4775 

12740 

5.324 

-2771 

-2668 

-2541 

-2619 

-2519 

-2393 

-1045 

4024 

11155 

9.500 

-2265 

-1872 

-1455 

-2130 

-1736 

-1317 

-  299 

3157 

7605 

16.500 

-  713 

-  587 

-  445 

-  747 

-  621 

-  479 

-  84 

1107 

2505 

Table  3.3c  Plate  beam  CX4 


Y 

oy  (psi) 

(in.) 

A1 

A2 

A3 

B1 

B2 

B3 

Cl 

C2 

C3 

0.676 

-48270 

-37590 

-26240 

-48610 

-37800 

-26415 

-49440 

-38595 

-27255 

3.000 

-43650 

-34230 

-24500 

-43315 

-33880 

-24220 

-43200 

-24330 

5.324 

-39980 

-30840 

-21740 

-40005 

-30960 

-22030 

-39965 

-31100 

-22435 

9.500 

-21720 

-21230 

-15420 

-27150 

-21320 

-15550 

-27130 

-21340 

16.500 

-  9105 

-  7235 

-  5325 

-  9110 

-  7207 

-  5277 

-  9160 

-  7225 

-  5270 

D1 

D2 

D3 

El 

E2 

E3 

FI 

F2 

F3 

0.676 

-50570 

-38710 

-27930 

-51640 

-39780 

-29005 

-12830 

+28290 

70310 

3.000 

-43655 

-33720 

-24340 

-44135 

-34210 

-24830 

-11605 

25560 

63520 

5.324 

-39870 

-31220 

-22580 

-39765 

-31120 

-22490 

-10430 

22850 

56830 

9.500 

-26980 

-21315 

-15520 

-26860 

-21195 

-15400 

-  7100 

16000 

38930 

16.500 

-  9136 

-  7210 

-  5315 

-  9162 

-  7236 

-  5343 

-  2400 

5389 

13220 
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Table  3.3d  Plate  beam  CX  5 


m 

Oy  (psi> 

A1 

A2 

A3 

B1 

B2 

B3 

Cl 

C2 

C3 

0.676 

-18990 

-11890 

-4678 

Bn 

-12480 

-4912 

-21430 

-13650 

-  5854 

3.000 

-18160 

-11120 

-4005 

-11220 

-4189 

-18680 

-11820 

-  4985 

5.324 

-10110 

-2875 

ns 

-10300 

-3573 

-17060 

-10880 

-  4788 

9.500 

-  7214 

-3338 

-11150 

-  7246 

-3347 

-11165 

-  7225 

-  3284 

16.500 

-  3494 

-  2229 

-  963 

-  3651 

-  2371 

-1094 

-  3780 

-  2493 

-  1212 

D1 

D2 

D3 

El 

E2 

E3 

FI 

F2 

F3 

0.676 

-23300 

-14550 

-6405 

-24480 

-15740 

-7608 

3672 

31575 

60930 

3.000 

-18720 

-11970 

-5518 

-19255 

-12515 

-6070 

3264 

28630 

55050 

5.324 

-16770 

-11000 

-5217 

-16655 

-10895 

-5118 

3011 

25750 

49160 

9.500 

-11200 

-  7256 

-3243 

-11070 

-  7123 

-3108 

2466 

18135 

33610 

16.500 

-  3653 

-  2411 

-1162 

-  3687 

-  2445 

-1196 

824 

6105 

11365 

Table  3.3e  Plate  beam  CX6 


Y 

0y  (psi) 

(in.) 

A1 

A2 

A3 

B1 

B2 

B3 

Cl 

C2 

C3 

0.676 

-11410 

-5070 

1312 

n 

-  6115 

461 

-14560 

-  7812 

-  1070 

3.000 

-11550 

-4975 

1625 

B 

-  5470 

906 

mm 

-  6619 

-  512 

5.324 

-11580 

-4742 

2110 

B9 

-  5015 

1183 

Bn 

-  5965 

-  464 

9.500 

-  6985 

-3453 

77 

-  7115 

-  3641 

-  167 

-  7194 

-  3776 

-  354 

16.500 

-  1936 

-  781 

371 

-  2351 

-  1243 

-  139 

-  2517 

-  1444 

-  379 

D1 

D2 

D3 

El 

E2 

E3 

FI 

F2 

F3 

0.676 

-17150 

-9085 

-1510 

-18270 

-10320 

-2770 

7100 

31300 

57310 

3.000 

-13360 

-7155 

-1157 

-13900 

-  7720 

-1740 

6500 

28450 

51670 

5.324 

-11400 

-6279 

-1159 

-11280 

-  6173 

-1070 

6020 

25725 

46180 

9.500 

-  7185 

-3829 

-  442 

-  7055 

-  3693 

-  301 

4514 

18140 

31640 

16.500 

-  2322 

-1220 

-  86 

-  2364 

-  1263 

-  130 

1573 

6189 

10650 
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CHAPTER  4 


SPECIAL  CLASS  OF  STRUCTURAL  PROBLEMS 
4.1  Introduction  to  Propeller  Blades 

Marine  propeller  blades  present  a  special  class  of  structural  problems  for  which  no 
generally  satisfactory  solution  is  presently  available. 

Screw  propellers  have  been  and  remain  the  principal  device  used  to  move  a  ship. 
Despite  their  importance,  there  is  little  or  no  realistic  approach  by  means  of  which  an  accu¬ 
rate  evaluation  of  propeller  stresses  can  be  obtained.  The  problem  lies  in  the  difficulty  of 
describing  a  blade  design  in  simple  mathematical  terms.  Until  recently,  existing  methods  ap¬ 
plicable  to  screw  propellers  have  relied  heavily  on  practical  experience  and  semi  theoretical 
considerations.  They  provide  a  criterion  of  stress  rather  than  the  actual  surface  stresses. 

In  the  past,  analytical  methods  for  predicting  blade  stresses  have  been  developed  by 
using  “beam”  theory  or  “shell”  theory.  The  use  of  elementary  beam  theory  was  first  pro¬ 
posed  by  Taylor  (1933)  who  treated  a  blade  as  a  cantilever  attached  to  the  propeller  hub  and 
recommended  that  stresses  be  calculated  for  cylindrical  blade  sections  with  the  neutral  axis 
parallel  to  the  nose-to-tail  (pitch)  line*  of  the  expanded  section.  Cantiliver  beam  theories 
have  yielded  reasonable  estimates  of  stresses  at  certain  selected  points  of  relatively  straight 
and  narrow  blades.  Some  modified  forms  of  beam  theory  have  been  proposed  (Rosengh,  1937; 
Hancock,  1942;  Schoenherr,  1963;  and  McCarthy  and  Brock,  1969)  for  wide-bladed  propellers 
with  blade  w'idth  to  length  ratios  of  about  1.  The  shell  theory  approach  was  first  proposed  by 
Cohen  (1955)  who  treated  a  simplified  propeller  blade  model  as  a  helicoidal  shell  with  vari¬ 
able  thickness  and  infinite  width.  However,  where  this  approach  was  applied  to  the  problem 
of  a  shell  of  finite  width,  it  was  impossible  to  produce  a  solution  to  satisfy  the  boundary  con¬ 
ditions.  Later  studies  included  those  of  Conolly  (1961)  and  others  (General  Applied  Science 
Laboratory,  1963  and  Atkinson,  1968).  Shell-type  theories  that  incorporate  broad  assumptions 
do  not  appear  to  offer  tangible  improvement;  moreover,  they  are  rather  involved  for  routine  de¬ 
sign  purposes.  Analytical  methods  that  attempt  to  predict  blade  stresses  on  the  basis  of 
conventional  mechanics  have  not  been  eminently  successful. 

Considerable  effects  have  been  devoted  to  measuring  blade  strains  on  both  model  and 
prototype  propeller  blades  (Connolly,  1961;  Wereldsma,  1965;  McCarthy  and  Brock,  1969; 
and  Boswell,  1969).  In  certain  cases,  good  agreement  was  obtained  between  beam  theory 
and  measured  data.  However  because  of  the  large  number  of  factors  involved,  care  must  be 
taken  in  drawing  general  conclusions  from  limited  measurements. 

*See  Comstock  (1967)  for  propeller  terminology. 
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The  trend  in  shipbuilding  to  full  afterbodies  for  mammoth  tankers  and  bulk  carriers  and 
to  higher  speeds  for  modern  naval  vessels  has  been  accompanied  by  large  irregularities  and 
fluctuations  in  the  wakes  of  such  ships.  As  a  result,  the  propeller  experiences  increased  dy- 
namic  excitation  and  generates  severe  vibratory  forces  on  the  hull  and  propulsion  system. 

Propeller-induced  vibration  is  one  of  the  main  problems  associated  with  ship  propul¬ 
sion  employing  the  screw  propeller.  The  thrust  derived  from  blade-lift  force  is  unsteady 
when  the  blades  rotate  in  a  nonuniform  velocity  field  behind  the  ship.  The  interaction  of 
these  unsteady  forces  with  the  hull  and  appendages  causes  the  excitation  of  the  ship  by  the 
propellers.  Blade  skew,  high  blade  area  ratios  (i.e.,  wider  blades),  and  a  larger  number  of 
blades  per  shaft  have  all  been  tried  in  an  attempt  to  reduce  vibration.  These  innovations  of 
propeller  geometry  drastically  alter  blade  displacement  patterns  (Dhir  and  Sikora,  1971,  1972) 
and  render  the  standard  methods  (i.e.,  beam  theory)  invalid.  If  blade  design  is  to  have  a 
sound  and  rational  basis,  then  an  effective  analytical  method  is  clearly  required  so  that  suit¬ 
able  blade  strength  and  stiffness  can  be  determined  for  a  specified  ship-operation  task. 

O 

The  finite  element  procedure  outlined  in  Chapter  3  will  now  be  used  to  analyze  a 
screw  propeller  in  its  more  general  form,  that  of  a  highly  skewed  propeller.  The  computed 
results  will  then  be  compared  with  measured  displacements  and  stresses  derived  from  mea¬ 
sured  strain  under  steady  pressure  loading.  This  study  will  provide  the  basis  for  further  ex¬ 
tending  the  procedure  until  it  is  eventually  able  to  take  unsteady  stressing  due  to  dynamic 
effects  and  fatigue  behavior  into  account. 

4.2  The  Geometry  of  Skewed  Propellers 

When  properly  designed,  skewed  propeller  can  be  used  not  only  to  significantly  reduce 
propeller-induced  vibrations  but  also  to  improve  blade  performance  against  cavitation,  etc. 
(Cummings  et  al.,  1972). 

The  general  characteristics  of  screw  propellers  and  their  typical  design  process  can 
be  found  in  standard  sources  (Comstock,  1967).  A  propeller  is  usually  described  by  a  set  of 
architectural  drawings  that  show  various  views  of  the  blade  and  a  table  that  gives  detailed 
dimensions  of  blade  sections  including  pitch  angles,  camber  ordinates,  thickness  ordinates, 
and  chord  lengths.  These  are  useful  input  for  propeller  manufacturers.  Application  of  the 
finite  element  procedure,  however,  requires  analytical  specification  of  the  blade  geometry, 
i.e.,  the  definition  of  a  suitable  set  of  coordinate  systems. 

The  highly  skewed  propeller  is  a  recent  innovation  and  its  geometry  is  more  involved. 
The  skewed  blade  illustrated  conceptually  in  Fig.  4.1  is  shown  in  Fig.  4.2  as  a  projection  on 
the  XY  plane.  Two  right-handed  coordinate  systems  share  the  common  Z  (=  z)  axis  (coincid¬ 
ing  with  the  centerline  of  the  shaft).  The  common  Z-axis  is  taken  positive  downstream  from 
the  propeller.  One  of  the  systems  is  the  global  rectangular  coordinates  (A,  Y,  Z),  a  ship-fixed 
frame  of  reference  whose  origin  lies  on  the  shaft  centerline  at  a  convenient  longitudinal 
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COORDINATES  FOR  DEFLECTION  MEASUREMENTS 
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location.  The  other  is  the  local  rectangular  coordinates  (x,y,Z)  attached  to  the  rotating  pro¬ 
peller.  The  y-axis  for  this  system  passes  through  the  midchord  of  a  blade  root  section.  The 
y-axis  is  tangent  to  the  blade  reference  line  and  is  also  known  as  the  generator  line  of  the 
blade.  The  angular  displacement  of  the  generator  line  6  is  measured  in  the  AT»plane  from 
the  global  axis  OX,  positive  in  a  counterclockwise  direction.  The  skew  angle  0S  is  the  pro¬ 
jected  angle  of  the  blade  measured  from  the  generator  line  oy. 

Fig.  4.3  provides  a  pictorial  view  of  the  basic  parameters  that  will  provide  a  complete 
description  of  the  blade  geometry.  The  blade  section,  which  lies  on  a  cylindrical  surface  of 
radius  p  (usually  expressed  as  percentage  of  blade  tip  radius  R),  is  shown  in  Fig.  4.4  in  de¬ 
veloped  view  and  with  more  details.  The  blade  section  has  a  reference  line  (also  known  as 
the  geometric  pitch  line)  which  makes  an  angle  <f>  with  the  a;y-plane.  The  geometric  pitch 
line  usually  lies  along  the  nose-tail  line  (curved  line  CDE  in  Fig.  4.3).*  This  line  joins  the 
the  leading  point  and  the  trailing  point  of  the  blade  section  cut  by  the  cylinder  (p  =  constant). 

Another  important  baseline  is  the  blade  reference  line  which  is  defined  as  the  line 
through  an  approximate  “datum”  point  at  each  radius  (p  <  R).  This  datum  point  is  usually 
on  the  nose-tail  line  and  will  generally  be  the  midchord  point  (Fig.  4.4). 

On  the  cylindrical  surface  (p  =  constant),  a  coordinate  system  (£x,  f2,  p)  is  defined.  The 
variable  ^  is  measured  from  the  datum  point  along  the  nose-tail  line  positive  in  the  down¬ 
stream  direction  (i.e.  toward  trailing  edge).  £2  is  normal  to  the  fj-line  and  points  in  the  up¬ 
stream  direction  (Fig.  4.3).  The  system  is  left-handed  to  permit  blade  section  geometry  to  be 
specified  in  the  conventional  two-dimensional  orientation  such  that  positive  camber  is  in  the 
positive  (^-direction.  Finally,  the  blade-section  profile  ordinates  (see  Fig.  4.4)  Ey  =  ( Ec  + 

Et)  and  E L  =  (Ec  -  Et)  for  -0.5 c  <  <  0.5c  can  be  expressed  in  terms  of  the  camberline 

function  Ec  (£j,p)  and  thickness  function  E (  (£j,  p). 

“Rake”  and  “skew”  have  been  implemented  frequently  in  modern  propeller  designs. 

They  represent  certain  specialized  departures  from  orthodox  blade  geometry  in  the  fore-  and- 
aft  and  transverse  projected  views,  respectively.  Unfortunately  use  of  the  terms  has  different 
meanings  for  different  people.  Therefore  there  is  a  need  to  have  precise  definitions  of  skew 
and  rake  because  they  are  essential  to  a  valid  structural  analysis. 

Skew  is  defined  as  the  sucessive  “transverse  displacement”  of  the  blade  sections 
along  their  respective  pitch  helices.  A  skew  angle  is  a  projected  angle  in  the  (transverse) 
AT-plane  between  a  radial  line  (generator  line,  Fig.  4.2)  from  the  centerline  of  the  propeller 
hub  through  midchord  of  the  blade  root  section  and  a  radial  line  through  midchord  of  the  blade 
section  in  question.  Each  section  of  the  blade  has  a  different  skew  angle;  the  skew  angle 
6  of  the  tip  section  is  commonly  used  as  the  measure  of  the  skew  of  the  blade.  Rake  IK  is 


*Space  curve  CDE  is  also  known  as  pitch  helix  which  is  generated  as  the  blade  rotates; 


is  the  pitch  angle. 
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LINE  CDE  IS  THE  NOSE-TAIL  LINE  LOCATED  ON  A  CYLINDRICAL  SURFACE  OF  RADIUS  p 
(THE  BLADE  IS  CUT  BY  A  CYLINDRICAL  SURFACE  CONCENTRIC  TO  THE  PROPELLER  AXIS,  LOOKING  AFT) 


Figure  4.3  Local  blade  coordinate  system  (x,  y,  Z )  and  (p,  6,  Z) 
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Figure  4.4  Developed  view  of  cylindrical  blade  section 


defined  as  the  successive  longitudinal  “displacement”  forward  or  aft  of  the  blade  sections 
from  the  radial  reference  line  (Fig.  4.3).  These  definitions  of  rake  and  the  skew  allow  the 
datum  points  along  a  blade  reference  line  to  be  described  analytically  by  ( p ,  6(p),  $  ( p ))  with 
respect  to  the  local  cylindrical  coordinates  (p,  6,  z).  Hence  the  blade  section  profile  and 
complete  blade  surface  can  be  constructed. 

4.3  Experimental  Data 

Dhir  and  Sikora  (1971)  used  holographic  interferometry  to  predict  the  static  deflections 
of  a  highly  skewed  marine  propeller  blade  model.  The  aluminum  blade  model  (Fig.  4.5)  was 
loaded  on  one  side  by  a  specially  designed  pressure  chamber  under  controlled  air  pressure 
(Boswell,  1969  and  McCarthy  and  Brock,  1969).  Applied  pressure  ranged  from  0.032  to  0.192 
psi.  This  holographic  technique  uses  laser  as  the  light  source  and  is  capable  of  measuring 
the  displacement  components  of  complex  surfaces  with  a  resolution  of  about  10  p  in. 

The  three-dimensional  displacements  were  derived  from  fringe  patterns  formed  in  the 
holographic  interferogram  (Fig.  4.6).  The  interferometric  fringes  on  an  object  were  related  to 
optical  path  differences,  and  the  double  exposure  technique  was  used  in  its  formation  (in  this 
case,  air  pressure  was  applied  to  the  blade  after  the  first  exposure  was  made).  Dark  and  light 
fringes  represent  contours  of  constant  spatial  blade  displacements.  As  the  applied  pressure 
increases  on  the  propeller  blade,  the  fringe  pattern  remains  essentially  unchanged  but  the 
number  of  interference  lines  increases.  The  appropriate  displacement  components  can  be  ex¬ 
pressed  as  a  linear  combination  of  the  “number”  of  fringes.  This  number  is  counted  from  the 
(fixed)  hub  to  the  point  of  interest  (obtained  through  certain  designated  observation  points  on 
the  hologram  plate)  and  is  multiplied  by  the  wavelength  of  light  used  (A  =  5145  A;  see  Dhir  and 
Sikora  (1971,  1972)  for  a  detailed  description).  Refined  displacement  data  are  obtained  by 
applying  the  least-square  approximation  technique. 

The  fringe  pattern  of  the  hologram  indicates  that  the  contours  of  constant  displacement 
line  cross  the  circular  lines  which  correspond  to  blade  sections  of  constant  radii  at  various 
angles.  It  is  concluded  that  the  deflections  of  the  blade  do  not  follow  those  predicted  by  the 
conventional  method  based  on  “beam”  theory. 

Strain  measurements  on  a  highly  skewed  model  propeller  blade  were  first  made  by  Boswell 
(1969).  The  blade  model  of  2014  aluminum  was  similar  to  that  shown  in  Fig.  4.5  except  that  it 
had  a  skew  of  120  deg.  Strain-gage  rosettes  were  placed  along  cylindrical  blade  sections,  five 
each  at  the  30-,  50-,  and  70-percent  radii  and  three  at  the  90-percent  radius.  The  rosettes  were 
so  oriented  that  both  tangential  and  radial  strains  were  measured.  For  every  gage  located  on  the 
blade  face,  there  was  a  corresponding  gage  on  the  blade  back. 

The  blade  was  subjected  to  air  pressure  loading  by  using  the  technique  described 
previously.  The  strain  distribution  was  found  to  be  radically  different  from  that  pteviously 
measured  on  unskewed  blades.  The  high  principal  stresses  occurred  in  a  relatively  narrow 
band  which  extended  from  near  the  trailing  edge  at  the  blade  root  to  near  the  leading  edge 
at  the  90-percent  radius. 
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To  obtain  an  experimental  assessment  of  the  effect  of  skew  on  blade  stress  distribu¬ 
tions,  Cumming  et  al.  (1972)  conducted  further  tests  on  a  series  of  one-blade  versions  of  pro¬ 
pellers  with  various  degrees  of  skew.  These  propellers  were  made  of  the  same  material  and 
tested  under  the  same  loading  conditions.  Stress  data  for  a  blade  with  72-deg  skew  were  ob¬ 
tained  for  comparison  with  the  present  study. 

4.4  Finite  Element  Analyses 

The  three-dimensional  curved  elements  described  in  Chapter  3  are  used  here  to  predict 
static  deflections  and  stresses  for  a  highly  skewed  propeller  blade.  The  blade  model  is  sub¬ 
jected  to  a  uniformly  distributed  pressure  applied  on  its  back  surface,  and  its  root  section  is 
assumed  to  be  completely  fixed.  This  loading  is  not  the  anticipated  service  load  distribution 
but  was  selected  to  match  the  experimental  conditions.  Nonuniform  loading  represents  no 
problem  to  the  program  and  the  influence  of  load  distribution  on  stresses  can  be  investigated. 

The  propeller  blade  used  in  the  current  analysis  is  one  in  a  series  of  skewed,  research 
model  propellers  that  form  part  of  a  larger  study  undertaken  to  determine  the  blade  strength 
and  other  performance  characteristics  of  highly  skewed  propellers.  Each  propeller  consists 
of  five  equally  spaced  1/23-scale  propeller  blade  models  attached  to  a  cylindrical  hub  (Fig. 
4.7).  The  version  selected  for  this  study  (Fig.  4.5)  is  a  single-blade  model  of  2014  aluminum 
with  a  tip  radius  of  6  in.  and  a  skew  angle  of  72  deg. 

The  geometry  of  the  propeller  blade  has  been  shown  to  be  rather  involved  (Section  4.2). 
Application  of  the  finite  element  method  requires  precise  coordinate  definition  of  hundreds  of 
discrete  points  that  cover  the  complete  top  and  bottom  surfaces  of  the  blade.  A  sizable  effort 
is  therefore  involved  in  the  calculation  of  coordinate  data.  Existing  computer  programs,  e.g., 
Brockett  (1972),  and  other  data  generation  routines  are  designed  to  furnish  the  pertinent  geo¬ 
metric  description  of  a  propeller  blade.  Nevertheless  it  is  the  responsibility  of  the  user  to 
discern  the  validity  of  his  input. 

The  coordinate  data  for  this  analysis  were  obtained  from  the  geometric  output  generated 
by  a  numerical  machine*  from  input  of  blade  section  parameters.  As  seen  from  its  plotted  pro¬ 
jection  (Fig.  4.8),  the  coordinate  data  are  for  a  blade  situated  in  a  position  144  deg  ahead  of 
the  blade  shown  on  Fig.  4.2.  The  geometric  data  have  been  checked  visually  and  by  graphical 
means;  however,  rigorous  error  analysis  has  not  been  attempted. 

Discrete  point  locations  on  the  blade  surfaces  were  chosen  for  the  element  mesh  net¬ 
work.  Figure  4.9  shows  that  these  points  were  along  20-,  30-,  40-,  50-,  60-,  70-,  80-,  and 
90-percent  blade  radius,  etc.  They  are  those  points  for  which  blade  section  geometries  were 
given  and  for  which  displacement  as  well  as  strain  measurements  were  made.  It  is  only  natural 
that  these  points  are  adopted  in  the  formation  of  the  element  mesh. 


*A  numerically  controlled  machine  used  for  model  cutting. 
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COORDINATES  FOR  FINITE  ELEMENT  CALCULATION 


Figure  4.8  XY-  and  YZ-projections  of  a  highly  skewed  propeller  blade 
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Three  elements  were  used  to  span  the  chordwise  dimension,  and  four  elements  were 
used  to  cross  from  the  blade  root  to  its  tip.  Available  coordinate  data  and  element  shapes 
into  which  the  blade  was  discretized  were  taken  into  consideration.  This  arrangement 
allowed  a  direct  comparison  of  computed  and  measured  deflections  and  immediate  interpre¬ 
tation  of  the  stress  calculation. 

To  a  large  extent,  the  effectiveness  of  a  finite  element  analysis  depends  on  the  appro¬ 
priate  element  arrangement  or  mesh  pattern  which  is  used  to  represent  a  given  structure.  The 
maximum  number  of  elements  is  generally  dictated  by  considerations  of  accuracy  and  cost. 

The  accuracy  of  the  computed  results  increases  as  the  number  of  elements  increases,  pro¬ 
vided  the  resulting  elements  are  not  so  numerous  as  to  exceed  computer  size  limitation  or  the 
number  of  computations  are  so  large  that  round-off  errors  become  a  dominant  factor.  On  the 
other  hand,  as  the  number  of  elements  increases,  the  cost  of  input  preparation  and  computa¬ 
tion  processing  also  increases.  Obviously  a  compromise  is  required  for  these  two  competing 
considerations.  There  is  no  easy  way  to  attain  an  optimal  mesh;  nevertheless,  numerical  in¬ 
stability  or  distortion  of  element  performance  can  always  be  avoided  by  the  discreet  exercise 
of  good  engineering  judgment.  (Some  basic  rules  for  element  mesh  layout  were  described  in 
Section  3.2.4.) 

After  element  meshes  were  generated  and  nodes  labeled,  the  elements  were  numbered 
in  succession  such  that  the  maximum  front  width  (see  Section  3.4.2)  was  optimized.  In  the 
finite  element  method,  each  element  is  defined  by  a  series  of  nodes  to  which  it  is  connected 
and  its  orientation  (with  reference  to  local  axes)  is  determined  by  the  numbering  sequence  of 
the  nodes  defined  in  the  incidence  table.  The  first  two  numbered  nodes  serve  to  fix  the 
direction  of  element  axis-£.  The  element  axis-7?  is  taken  in  the  direction  of  a  line  from  the 
first  to  the  third  numbered  node.  The  £-axis  is  obviously  the  remainder  that  completes  the 
right-handed  element  axis  system.  In  the  case  of  propeller  blades,  nodes  can  be  placed  along 
a  cylindrical  surface  so  that  the  direction  of  local  element  axes  f  and  7?  will  respectively 
correspond  to  the  radial  and  tangential  axes  of  the  cylindrical  coordinates. 

Note  that  the  experimentally  determined  blade  displacements  were  measured  with  re¬ 
spect  to  the  global  coordinates  ( X ,  Y ,  Z)  when  the  blade  was  oriented  in  the  position  shown 
in  Fig.  4.2,  where  6 1  =  124.96  deg.  The  finite  elements  coordinates  obtained  from  the  numer¬ 
ical  machine  provided  a  blade  with  the  projection  shown  in  Fig.  4.9,  i.e.,  6  j  was  approximate¬ 
ly  269  deg.  Finite  element  calculations  gave  displacement  components  u  and  v  that  were  144 
deg  out  of  phase  (At^)  with  reference  to  measured  displacements;  accordingly,  the  element 
displacements  were  transformed  so  that  they  could  be  compared  on  the  same  scale  (Fig.  4.10). 

Computed  displacements  were  obtained  along  the  cylindrical  blade  section  at  30-,  40-, 
50-,  60-,  70-,  80-,  and  90-percent  radii  in  addition  to  points  near  the  blade  tip.  Measured  de¬ 
flections  by  holographic  interferometry  are  available  at  50-,  70-,  and  90-percent  radii  on  blade 
face;  these  are  plotted  on  Fig.  4.10. 
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Figure  4.10b  Displacement  V 


INCH  X  10~3  50  PERCENT  RADIUS 
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Computed  stresses  in  global  coordinates  at  interior  points  as  well  as  surface  points 
are  given  along  cylindrical  blade  sections.  Stress  outputs  are  available  at  30-,  50-,  70-,  and 
90-percent  radii  and  at  many  intermediate  points.  Stress  transformation  is  applied  so  that 
there  is  an  option  to  print  out  stresses  in  local  cylindrical  coordinates.  Stresses  in  most 
areas  are  obtainable  by  interpolation.  Computed  stresses  are  plotted  on  Fig.  4.11  together 
with  stresses  calculated  from  experimentally  obtained  strain  data. 

4.5  Discussion  of  Results 

Fig.  4.10  shows  that  blade  displacement  component  w  (parallel  to  the  centerline  of  the 
propeller  shaft)  was  generally  higher  near  the  trailing  edge  than  near  the  leading  edge  for 
various  blade  radii.  This  was  also  visible  from  the  interferometric  fringes  (Fig.  4.6)  repre¬ 
senting  displacement  contours.  This  indicates  that  displacements  for  a  skewed  propeller  are 
no  longer  parallel  to  the  cylindrical  blade  sections.  Clearly  this  is  a  case  for  which  the 
assumption  made  by  the  conventional  method  of  blade  design  is  not  valid.  Further  it  is  noted 
that  the  values  for  u  and  v  are  appreciably  larger  than  the  displacements  that  can  be  account¬ 
ed  for  by  beam  theory.  A  new  approach  is  needed  to  predict  the  elastic  behavior  of  a  skewed 
blade. 

Genalis  (1972)  treated  the  skewed  blade  as  a  thin  shell  and  idealized  it  as  an  assem¬ 
blage  of  discrete,  flat,  triangular  elements;  see  Fig.  4.12.  Because  of  the  complexity  of  blade 
profiles  and  blade  surface  curvatures  (such  as  described  previously  in  Figs.  4.3  and  4.7),  a 
large  number  (about  270)  of  these  shell  elements  is  required  in  order  to  achieve  a  reasonable 
approximation  to  the  geometry  of  a  skewed  blade.  Thus  a  substantial  effort  is  required  for 
input  preparations  including  all  coordinates  and  load  data.  (It  is  conceivable  here  that  the 
chance  for  input  error  is  increased  because  of  the  huge  amount  of  complex  input  data  involved.) 
Furthermore,  the  proportions  of  the  elements  in  the  region  near  the  hub  are  well  beyond  the 
limits  of  medium  thick  plate  or  thin  shell  theory  from  which  such  elements  are  derived.  Com¬ 
puted  displacements  were  given  for  the  shell  solution  along  cylindrical  sections  at  50-  and 
70-percent  radii  (Fig.  4.10).*  The  discrepancy  of  the  numerical  results  as  obtained  from  the 
flat  thin  shell  element  solution,  from  the  isoparametric  finite  element,  and  from  the  experi¬ 
mental  methods  can  be  attributed  (a)  to  the  geometric  mismatch  introduced  by  using  a  flat 
shell  element  approximation  for  this  highly  curved  blade  of  varying  thickness  and  (b)  to  the 
breakdown  in  the  suitability  of  thin  shell  theory  in  the  region  near  the  hub.  (Furthermore, 
some  inaccuracy  may  be  attributable  to  rounding  errors  developed  during  the  solution  process.) 
The  error  accumulation  evident  in  the  thin  shell  element  values  at  the  70-percent  radius  was 
even  more  pronounced  at  90-percent  radius. 


Any  inaccuracy  of  input  data  will  almost  inevitably  show  up  in  the  numerical  results.  The  reliability  of  com¬ 
puted  values  may  be  increased  if  an  appropriate  algorithm  is  available  to  screen  or  to  smooth  input.  The  toler¬ 
ance  in  the  fabrication  of  the  blade  model  was  said  to  be  within  5  percent 
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Figure  4.12  Flat  shell  element  representation  of  a  72-degree  skewed  propeller  blade 


The  changing  profiles  of  the  successive  sections  from  blade  root  to  tip  indicate  that 
the  skewed  blade  has  a  more  complicated  geometry  than  an  ordinary  shell.  The  blade  is 
treated  in  this  study  as  a  solid  continuum  of  general  shape.  It  is  represented  by  12  curved, 
three-dimensional  elements  which  are  capable  of  fitting  smoothly  in  any  complex  surface. 

Use  of  these  isoparametric  elements  gives  a  marked  improvement  in  the  numerical  results 
(Fig.  4.10).  Note  that  the  computed  displacements  agree  well  with  holographic  measurements, 
especially  in  displacement  component  w,  which  is  by  far  the  dominant  factor  of  the  total  blade 
displacement. 

The  largest  disagreement  between  the  displacements,  which  are  solved  for  directly  by 
the  two  procedures  (experimental  and  isoparametric  finite  element)  is  in  the  horizontal  dis¬ 
placement  normal  to  the  shaft  axis.  Possible  causes  for  such  a  discrepancy  might  be  slippage 
of  the  test  specimen  at  the  support  bolts  in  the  hub  and  mismatch  of  the  isoparametric  propel¬ 
ler  surface  with  that  of  the  test  specimen.  However,  since  the  differences  do  not  show  up  in 
the  other  displacements,  the  latter  explanation  seems  unlikely.  Furthermore,  support  slippage 
(i.e.,  a  difference  between  the  ideal  clamped  condition  of  the  finite  element  method  and  the 
fixity  achieved  in  the  test)  also  seems  unlikely  because  of  the  low  pressures  involved.  If 
movement  did  occur,  it  would  have  to  be  rigid-body  movement  normal  to  the  axis  of  the  shaft 
rather  than  the  more  likely  rotation  or  tipping.  The  final  possibilities  are  errors  in  reading 
and/or  possible  binding  of  the  blade  model  on  the  pressure  seal. 

The  stress  distribution  of  highly  skewed  propellers  is  expected  to  be  unsymmetric  and 
different  from  that  of  an  unskewed  blade.  The  maximum  radial  and  tangential  stresses  are 
located  away  from  midchord  and  toward  the  trailing  edge  in  the  main  body  of  the  skewed  blade. 
The  magnitude  and  distribution  of  stresses  along  the  cylindrical  chords  of  the  blade  varies 
with  blade  skew.  Computed  stresses  were  obtained  from  the  same  element  mesh  (Fig.  4.8) 
and  were  shown  together  with  experimental  stresses  in  Fig.  4.11.  No  calculated  stress  was 
given  by  the  shell  element  program  (Genalis,  1972). 

A  finer  mesh  would  have  given  a  closer  approximation  to  the  true  stress  distribution 
since  convergence  is  assured  with  the  compatible  elements  used  in  the  current  study. 

In  some  of  the  experiments,  difficulties  were  encountered  in  providing  a  precise  fit  be¬ 
tween  the  test  pressure  chamber  and  the  propeller  blade;  thus  free  movement  of  the  blade  edge 
was  not  always  obtained.  (It  was  also  noticed  in  the  course  of  stress  evaluation  for  a  series 
of  test  propellers  with  varying  degree  of  skews  that  maximum  stresses  and  blade  skews  did  not 
increase  in  direct  proportion.)  It  is  desirable  therefore  to  implement  a  comprehensive  test  pro¬ 
gram  to  enhance  the  reliability  of  experimental  results.  The  stress  agreement  shown  in  Fig. 
4.11  was  much  better  on  the  back  or  loaded  face  than  on  the  front  face,  and  the  agreement  was 
better  for  tangential  than  for  the  radial  stress.  However,  since  the  stresses  are  one  step  re¬ 
moved  from  direct  measurement  the  agreement  is  not  out  of  line  even  with  these  differences. 
Since  the  near  tip  region  of  the  blade  is  a  sensitive  region,  a  higher  density  of  elements  in  that 
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region  seems  desirable.  Coupled  with  the  comprehensive  test  program,  this  should  assist  in 
resolving  the  differences  and  provide  an  accurate  delineation  of  the  stresses. 

The  present  procedure  can  generate  an  estimated  complete  stress  field  for  the  propel¬ 
ler  and  thus  provide  a  detailed  account  of  the  stress  distribution  throughout  the  elastic  body. 
This  is  of  great  practical  value  for  a  structure  of  unconventional  configuration  where  conven¬ 
tional  methods  offer  little  guidance  as  to  the  pattern  of  stress  distribution  for  a  structure 
like  a  skewed  propeller.  A  development  that  utilizes  three-dimensional  curved  elements  is 
probably  the  most  effective  analytical  method  currently  available  and  is  capable  of  represent¬ 
ing  both  the  complex  blade  geometry  and  load  distributions  with  ease  and  validity. 
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CHAPTER  5 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  rational  approach  established  herein  represents  the  first  step  in  a  general  solution 
of  the  propeller  strength  problem.  A  finite  element  displacement  model  is  utilized  to  predict 
the  behavior  of  an  elastic  body  of  arbitrary  shape  under  static  loads.  Compatible  solid  ele¬ 
ments  in  their  general  form  are  adopted.  This  formulation  bypasses  the  constraints  of  simpli¬ 
fying  assumptions  and  allows  a  closer  approximation  to  the  true  structural  configuration  than 
is  possible  by  most  other  approaches,  for  instance,  by  classical  plate  or  shell  theories.  Solu¬ 
tions  are  subsequently  obtained  for  displacements  and  internal  stresses. 

The  numerical  results  obtained  by  the  refined  curved  elements  described  in  Chapter  3 
are  distinctly  superior  to  those  obtainable  with  commonly  available  simple  triangular  or  rec¬ 
tangular  finite  elements.  The  displacements  obtained  for  a  complex  structure  under  prescribed 
loading  had  a  good  degree  of  precision  even  with  the  relatively  coarse  mesh  employed.  This 
study  confirms  the  findings  of  Clough  (1969),  Zienkiewicz  (1971),  etc.  that  significant  improve' 
ments  in  element  performance  are  obtainable  when  higher  order  displacement  functions  are  em¬ 
ployed  in  the  element  formulation. 

For  structures  loaded  so  that  a  steep  stress  gradient  prevails,  stress  predictions  can 
always  be  improved  by  using  progressingly  finer  element  meshes.  The  element  presented  in 
this  study  showed  a  rapid  rate  of  convergence  as  the  mesh  was  refined. 

The  mapping  technique,  which  is  based  on  the  application  of  the  shape  function  as  de¬ 
scribed  in  Chapter  3,  provides  an  expedient  way  to  describe  the  geometry  of  a  general  continu¬ 
ous  surface.  Load  computations  for  a  prescribed  pressure  distribution  are  therefore  conducted 
in  a  routine  manner. 

The  method  developed  in  this  study  is  particularly  valuable  to  structures  with  arbitrary 
configurations  such  as  skewed  or  unsymmetric  geometry.  Large  savings  in  labor  and  comput¬ 
ing  costs  are  possible  when  the  flat-plate  or  shell-element  representation  of  a  complex  skewed 
body  is  replaced  by  curved  elements.  The  present  development  offers  perhaps  the  only  real¬ 
istic  solution  available  for  many  difficult  structural  problems  of  a  three-dimensional  nature, 
such  as  a  highly  skewed  propeller.  The  method  now  permits  rapid  solutions  for  structures 
with  varying  degrees  of  skew,  camber,  and  loading. 

Among  other  significant  applications,  the  present  procedure  can  provide  stress  calcula¬ 
tions  at  the  root  section  of  a  propeller  blade  joining  the  hub,  at  the  complicated  intersection 
of  cylindrical  walls  such  as  a  piping  joint,  and  at  the  delicate  interface  area  of  a  sandwich 
construction.  Further  extension  of  the  analysis  may  include  considerations  of  anisotropic 
material  properties  as  well  as  the  nonlinear  or  plastic  behavior  of  materials  (Gupta,  1971). 

The  current  method  provides  a  general  and  realistic  solution  to  the  structural  problem 
of  a  propeller  under  static  loads,  and  eventually  a  procedure  can  be  evolved  for  the  complete 
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solution  of  propeller  strength  under  all  possible  service  loads.  Inasmuch  as  a  propeller 
blade  experiences  varying  vibratory  forces  in  addition  to  the  “nearly  steady”  thrust  that  is 
always  present  during  ship  operation,  a  dynamic  analysis  that  includes  an  investigation  of 
fatigue  behavior  of  the  complex  ship  component  is  naturally  of  prime  significance. 

It  is  recommended  that  the  current  procedure  be  extended  to  incorporate  a  mass  matrix 
of  the  propeller  blade  and  an  auxiliary  matrix  to  account  for  the  effect  of  “added  mass”  of 
the  fluid  medium  surrounding  the  blade.  When  a  pressure-time  relation  for  a  propeller  is  pre¬ 
scribed,  a  complete  dynamic  response  can  be  obtained  by  using  a  time  increment  numerical 
technique.  The  effect  of  damping  may  also  be  incorporated. 

A  sound  propeller  design  is  extremely  important  to  the  successful  operation  of  a  ship 
and  has  a  significant  bearing  on  its  performance.  It  is  worth  noting  that  improved  design 
methods  including  a  dynamic  analysis  are  being  vigorously  pursued  by  Norwegian  and 
Japanese  ship  builders.  To  remain  competitive,  U.S.  ship  builders  must  take  advantage  of 
new  technologies  and  be  instrumental  in  their  development. 
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APPENDIX 

EXAMPLE  OF  FORTRAN  PROGRAM  FOR  NUMERICAL  CALCULATION 
OF  AN  ELEMENT  STIFFNESS  MATRIX 

An  outline  of  the  process  is  described  in  Section  3.2.2. 
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c>o  ooooooooooooooo 


LIST 

SUBROUTINE  STIFFB(XYZ,E,GNU,cK) 

COMMON  /  V  A  D  7 /  I NTGl R , ND , NPT , NPT S , NPZ 
COMMON  /BLAV/GM( 3, 5) ,PX(3> ,PY(3> ,0 
SUBROUTINE  TO  CALCULATE  THE  STIFFNESS  MATRIX  SK  FOR  ISO-PARAMETRIC 


S  K= INTE  GRAL  PHI  DV 
WHERE  PHI=  B  T  X  0  X  B 


X.  .  . 
0.  . . 

0  •  •  » 

Y  .  .  . 

0 «  . . 

0.  .  . 

X . . . =VcCT  OR  OF  PARTIALS 
WITh  RESPECT  TO  X 

OF 

N  (I) 

0- 

u. .  . 

Y.  .  . 

1)  .  «  b 

X  •  •  » 

Z.  .. 

0  •  *  • 

Y...=VECTOR  OF  PARTIALS 
WITH  RESPECT  TO  Y 

OF 

N  ( I ) 

J  •  .  . 
Z.  .  . 

z. .  . 

0  •  •  • 

Y.  .  . 

X.  .  . 

Z...=VECTOR  OF  PARTIALS 
WITH  RESPECT  TO  Z 
o. . .  =  :  VECTOR 

OF 

N  ( I ) 

NOTE 

PRIOR 

THAT 

TO 

THIS 

CODING 

MATRIX 

. 

MULTIPLICATION  HAS  BEEN  DON 

ALGEBRAICALLY 

REAL  CHI  (20  *3)  ,  CHI  X  ( 2  U  )  ,  CH I Y  (2 1. )  ,CHIZ(20) 

EQUIVALENCE  ( CHI  X < 1 )  ,  CH  I  < i  ,  1 )  )  ,  { C H  I Y  ( 1 >  ,  CH I ( 1 , 2 ) )  ,  (  CH IZ ( 1 ) , CH I (1 , 3 
X)  ) 

REAL  JAC  (3 ,3)  ,  DcT  J  ,  Jli  ,  J12  ,  Jl3  ,  J2i  ,  J22  ,  J23  ,  J31 ,  J32  ,  J33 
EQUIVALENCE  (JAC(1,1),J11),(JAC(1,2),J12>,(JAC<1,3>,J13>,(JAC<2,1> 
X,J21)  ,  (JAC (2,2), J2  2)  ,  (JAC (2, 3)  ,JZ3)  ,  (JAC (3,1)  , J31) ,  (JAC (3, 2)  ,J32)  , 
X  ( JAC (3  ,  3)  ,  J33) 

DI  Mt  NS  ION  E  K  (  6  U  ,  b  0 )  ,  SK(6  (l,6u)  , XYZ  (  20 ,3) 

DIMENSION  A (36) , H(36)  ,ISWK(6u>  ,TZ (3,3)  >S  (13) 

EQUIVALENCES  ABOVE  ARE  LOGICALLY  NECESSARY 
EQUIVALENCES  BELOW  ARE  USED  ONLY  TO  CONSERVE  SPACE 
EQUIVALENCE  ( S (1 ) , JAC (1 , 1) )  ,  (  S (1 )  , SI ) , ( S ( 2 ) , S2)  , ( S ( 3)  ,S3)  ,  ( S( 4) , S  1 
XXXF'l)  ,  (S(5)  ,SlYYPl)  ,  (S(6>  ,S2YYPl)  ,  <S  (7)  ,S2ZZPi>  ,  <S  (7)  ,S3ZZP1)  ,  <S  <  0 
X)  ,ZZ)  ,  (S (1 0 > ,CHIXI >  ,  (S(l 1) ,CHIYI >  , (S (12) ,CHIZI)  , (S ( 13) ,0ETJ) ,  <S( 10 
X)  , SUM)  ,  <S( 1U>  ,CON> 

C 

C  GAUSSIAN  QUADRATURE  CONSTANTS 

DATA  A/ 


1 

U  .  0  0  J  l)  II 

o  ouou 

U  DUO, 

-0. 57735 

02691 

3963, 

0.57735 

02691 

0963  , 

3 

-  IJ  •  7  7  4  5  9 

66b92 

4  14  0, 

G  .  U  6  U  U  U 

Jo  U  O  L 

0  0  0  j  , 

C . 77459 

66692 

4148  , 

4 

-0.86113 

63115 

9  4  0  5, 

-0 . 33998 

1  j  4  35 

8436  , 

0.33990 

10435 

04  36, 

4 

0.86113 

u  3  1 1  5 

9  4.5, 

-C . 9. bl 7 

98  4  5  9 

30  66  , 

-  0 «  5  30  4  6 

93141 

0  5  60  , 

5 

U  .  0  0  0  U  0 

OUUOO 

0  OU  U, 

0 • 53846 

93101 

0566, 

0.90617 

93459 

3  8  66, 

6 

-0.93246 

9  5 1 4  2 

V  315, 

- u . 66120 

93864 

6  b  2  7  , 

-0.23651 

91860 

8  3  20  , 

6 

J  .  2  3 :)  b  1 

9  1  3  o  j 

0  323  , 

u . 66120 

93  8  64 

6627  , 

0.9324b 

95142 

0  315, 

7 

-0.9491G 

79123 

4  27  o  , 

-C . 74153 

11855 

9939  , 

-0. 4U58  4 

51613 

7740  , 

7 

U  .  0  0  0  0  0 

u  o  u  o  u 

ll  0U  0  , 

0 . 4  U  5  8  4 

51513 

7  7  4  u  , 

0.74153 

11855 

9939  , 

7 

0  .949  1  0 

79123 

4  276, 

-C  .  9 6  u  2  3 

9o5  64 

9  7  54  , 

-o. 7  9  o  6  fa 

647  74 

13  63  , 

0 

-  0 .92553 

24099 

1  633  , 

-  0  .  10  34  3 

4&4  2 4 

9  5  65, 

0 . 13343 

4  6  •■+2  4 

9  5  65, 

0 

3.525  53 

2  4099 

1  63  3  , 

0  .  7  9  b  o  6 

047  74 

1363  , 

U  .  9  6  0  2  8 

98  564 

9754/ 

DATA  H/ 

1 

2  .  C  C  5  u  C 

C  3  J  0  0 

’J  0  o  J  f 

1  .  (  l.  J  U  J 

J  U  0  u  C 

j  0  J  b  j 

X  •  O  J  O  J  0 

0  U  0  j  0 

u  0  0  , 

3 

J  .  5  5  5  5  5 

5  5555 

5  556, 

1/  •  tC  0  y  0 

68886 

6  8  3  9  , 

o .  5  55  5  5 

5  5  555 

5556  , 

4 

3 .34785 

4  84  5  1 

3  7-,  5  , 

.C  .  fc 5 2 1  h 

615  4  0 

6  255  , 

U  .  6  5  2  i  4 

516  4  8 

62  55, 

4 

U  .347  8  5 

40451 

3  745, 

0. 23o92 

68850 

5619, 

U. h7&62 

66704 

9937, 

5 

3  •  6  6  8  8  i  > 

8  0888 

8  889  , 

C  •  m  7  8  -o  2 

86  7  4 

9937, 

0.23692 

b  8  8  5  j 

5619, 

6 

J.  17132 

44923 

7  917, 

3  .  3  6  v.  7  6 

1 5  7  3  C 

4014, 

J  a  4  6  79  1 

39345 

7269  , 

6 

0 . 4  o  7  9 1 

39345 

7  269  , 

C  .  3  6  j  7  6 

157  3  0 

4014  , 

U. 17132 

4h92  3 

7917, 

7 

0.12940 

49661 

6  8C7, 

0. 2797U 

539  14 

3923  , 

0 . 38103 

0  0  50  5 

0512  , 

7 

C  .  4  1  7  9  5 

91836 

73  4  7  , 

C  .38163 

0  0  5  o  5 

0512, 

o.  2797  0 

53914 

8920  , 

7 

3.12940 

49bol 

o  867, 

0 .  IE  122 

653  62 

9  0  38  , 

0.22238 

10  344 

5337, 

8 

0 . 3 1 3  7 1. 

6  0458 

7  709  , 

E  .  3  fc  2  6  8 

370  33 

78  36  , 

0.36268 

37033 

7836, 

0 

U.  313  70 

b  &4  5  8 

7  70  9, 

0. 22233 

103  44 

5337  , 

0.10122 

05362 

9038/ 

C  SIZE  CONSTANTS 
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U  f'l  A  NtJr^/2-J/,LNPTS'/4u/,KSIZt/6J/>KSX^tM/59/ 

C  CHECK  .MPT  A  NO  COP  PEC  T  IF  NECESSARY 
IF  <  NPT . GT. 3 ) N  P  T  = 3 
IF (NPT.LT.1)NPT=1 
IPT  =  NP  T*(NPT-l)/2 
IPZ  =  NPZ*  (  NPZ-l)  /2 
C  ZERO  MATRIX 

CALL  ZLFOZ (GM, 10 > 

CALL  ZtROZ <SK,3b JC> 

AA-l.-GNU 

m=gnu 

CC  = ( i. -2. ’GNU) /2 . 

DK  =  E  /  (  (  1.  4- GNU)  *(  1.-2.  *GNU>  ) 

NPZi=  NPZ  +  1 
N  D  Z  2  =  NPZ+2 

C  LOOFS  OVER  XX, YY,  AND  ZZ  FUR  GAUSSIAN  QUADRATURE  OVER  VOLUME 
00  2.  I XX= 1 , NPT 
XX-A (IXX+IPT) 

HXX=H( IXX+IPT) 

DO  23  I YY  =  1 , NPT 
YY-A (I YY+IPT) 

H X  X H  Y  Y  =  H  < X *H  (  IYY+IPT) 

00  22  I Z Z= 1 , NPZ2 
11-  A( IZZ+IPZ) 

IF(IZZ.EQ.NPZl)  ZZ=-1.G 
IF (IZZ.EQ.NPZ2)  ZZ=1.0 

C  CALCULATE  pARTIALS  OF  N 1 1 )  WITH  RESPECT  TO  XX,  YY,  AND  11 
DO  32  1=1,2 
S1=2*I -D 
S1XXP1  =  S1*XX+  1. 

S1YYP1=S1*YY+1. 

00  33  J=l, 2 
S2=2»J-o 
JP=2*C J-l) +1 
32YY?i=S2¥YY+l. 

S2ZZPi=S2¥ZZ+l. 

CHlX<JP  +  i3)=31¥U.23*S2YYPi*(l.-ZZ**2) 
CH1Y(JP+5)=S2¥3.23¥S1XXP1¥<1.-ZZ¥¥2> 

CHIZ ( JP  +  rt>  =  -3.5¥ZZ*SlXXPl»S2YYPl 
CHIX(JP+l2)=Sl*j.25*S2ZZPl*(l.-YY*»2) 

CHIY(JP+12) =-n.5*YY*SlXXFi*S2ZZPl 
CHIZ (J P+12) =3 2*3. 2 5* S1XXP1¥(1.-YY¥¥2> 

CHIX  (JP  +  16)  =-0.‘3»XX»SlYYPi»S2ZZPl 
CH1Y ( JP  +  lt>) =Si*o .25*S2ZZP1*(1.-XX**2) 

CHIZ  ( J  P  +  lc ) =  S  2  ¥U .25*S1YYP1*  ( 1 . -X  X  ¥¥2 ) 

DO  3  3  K  = 1 , 2 

S3=2*<-3 

KP  =  4  ¥ ( K-l)  +JP 

S3ZZP1=S3¥ZZ+1. 

CHIX(KP)=S1¥3.125¥S2YYP1¥S3ZZP1¥(S1¥Z,¥XX+S2¥YY+S3¥ZZ-1.) 
CH1Y (<P) =S2¥3 . 12  5¥S1XXP1¥S3ZZP1* (S2*2. *Y Y +S 1¥ XX +S 3* 11- 1 . ) 

33  CHIZ  UP) =  S3¥0 . 12 :>* SI XXP1 ¥S2 Y Y PI* (S3*2.*ZZ+S1*XX+S2*YY-1.) 

C  CALCULATE  JAC=PARTIAlS  X  XYZ 
00  4L  1=1,3 
DO  ■+  *  J=  1 ,  o 
SUM=U. L 

00  41  K=1 , NPTS 

41  S JM=  SU M  +CH I  (K,I)  *XYZ(K,J) 

40  JAC ( I , J) =SUM 
DO  42  L  =  1 , 3 
PX ( L ) =  JAC ( 1 , L ) 

42  PY ( L ) =  JAC ( 2 , L ) 

CALL  SURFN(TZ) 

C  INVERT  JAC  AND  CALCULATE  DET  JAC 
CALL  MINVDPtJAC, 3,3,QETJ,I0, IEXP) 

IF(ID-l)  31,31,33 
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33  PK1N  f  91,  ID 

91  FORMAT  <1X,I1Q,*TH  ORDER  PRONCIPAL  MINOR  IS  ZERO*  > 

C 

STOP 

31  CONTINUE 

C  CALCULATE  PATIALS  WITH  RESPECT  TO  X,  Y,  AND  Z  FROM  PARTIALS  WITH 
C  RESPECT  TO  XX,  YY ,  AND  ZZ  BY 
C  NEW. PARTI ALS=JAC  INVERSE  X  OLD. PARTIALS 
DO  50  I =1 , NPTS 
CHIX I=CHIX (I) 

CHIYI=CHIY (I) 

CHIZI=CHIZ (I) 

CHIX  (I ) =J11*CHIX 1+ J12*CHIYI+J13*CHIZI 
CHIY (I >=J21*CHIX  I+J22*CHIYI+J23*CHIZI 
50  CHIZ (I ) =J31*CHIX I+J32*CHIYI+J33*CHIZI 

WRIT E ( 4 )  IXX, IYY  ,IZZ , (CHIX (I)  ,CHIY (I)  ,CHIZ  CI>  ,1  =  1 , NPTS) ,11,11 
IF  (IZZ.GT.NPZ)  GO  TO  22 
C 

C  CALCULATE  INTEGRATION  CONSTANT 

CON=OETJ*DK*HXXHYY*H(IZZ+IPT) 

C  DO  SUMMATION  FOR  HALF  OF  MATRIX 

00  60  1=1, NPTS 
X1=CHI X  ( I) 

Y1  =  CHI Y  (I) 

Z1=CHIZ (I) 

DO  60  J= 1 , NPTS 
X2=CHI X ( J) 

Y2=CHI Y  (J) 

Z2  =  CHIZ  ( J) 

SK<I ,J+NPTS)=SK( I, J+NPTS)+C0N*<3B*X1*Y2+CC*Y1*X2) 

SKt I, J+LNPTS)  =SK  (I , J+LNPTS) +CON*  (  86*X1*Z2  *-CC*  Z1*X  2) 

SK (I+NPTS, J+LNPTS) =SK(I+ NPTS, J+LNPTS) +  CON* (BB*Y1*Z2+CC*Z1*Y2) 

IFU.LT. I)  GO  TO  61 

X2  =  X 1*  X  2 

Y2=Y1* Y2 

Z2=Z1*Z2 

SK<I , J) =SK(I, J)+CON* < AA*X2+CC*(Y2+Z2> ) 

SK (I+NPTS, J+NPTS )=SK (I+NPTS , J+NPTS) + CON* ( AA* Y2+ CC* ( X2 +Z2 ) ) 

SK (I +L NPTS .J+LNPTS )  =  SK( I +LNPTS, J  +  LNPTS) +CON* ( A A*Z 2  +  CC* ( X 2  +  Y 2)  ) 

61  CONTINUE 

60  CONTINUE 

22  CONTINUE 

CALL  GLOAD  ( X X ,Y Y , HX XHY Y ,X Y Z) 

20  CONTINUE 

RETURN 
END 
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